A&A 366, 945-964 (2001)
DOI: 10.1051/0004-6361:20000353

Falling Evaporating Bodies around Herbig stars. A theoretical study

H. Beust - C. Karmann - A.-M. Lagrange

Laboratoire d'Astrophysique de l'Observatoire de Grenoble, Université J. Fourier, BP 53,
38041 Grenoble Cedex 9, France

Received 19 April 2000 / Accepted 24 November 2000

Abstract
Transient spectral absorption events monitored now for years towards the star $\beta\:$Pictoris have been interpreted as resulting from the transit across the line of sight of evaporating star-grazing kilometer-sized bodies (Falling Evaporating Bodies, or FEBs). Several Herbig Ae/Be stars of various ages have been observed to exhibit somehow similar absorption events that have been attributed to similar FEB events. We investigate here this question from a modeling point of view. Adapting the FEB simulation code we had developed earlier specifically for $\beta\:$Pic to the case of typical Herbig Ae/Be stars, we try to derive in which conditions FEB-like objects may generate detectable transient absorption events. We compare these conditions with those found in the case of $\beta\:$Pic. A major difference with $\beta\:$Pic is that Herbig Ae/Be stars have strong stellar winds (10^-9- $10^{-7}~M_\odot$yr^-1). Those winds appear to have a drastic interaction with the gaseous material escaped from the FEBs. With the presence of such stellar winds, the spectral signatures of FEBs are not detectable, unless their mass loss rate is huge. This translates into very large bodies ($\sim$100 km size), instead of $\sim$15 km for $\beta\:$Pic FEBs. This appears unrealistic in terms of amount of planetesimal mass needed in the disks surrounding these stars. We discuss then the validity of the FEB hypothesis for specific example stars. It turns out that for the younger (a few $10^6\,$yr old) Herbig Ae/Be stars like AB Aur, with well identified winds yr^-1, the variable features sometimes observed are not likely to be due to FEBs, unless produced in wind free cavities. For older ( yr old) stars Herbig Ae/Be like HD 100546, the FEB scenario could still explain the spectral events observed, but either the wind must to be weaker than $\sim$10 $^{-10}~M_\odot$yr^-1 (which cannot be excluded so far), or the FEBs approach the star in wind free cavities.

Key words: stars: circumstellar matter - stars: early-type - stars: $\beta\:$Pic - methods: numerical -
techniques: spectroscopic

1 Introduction. The Falling Evaporating Bodies scenario among various stars

1.1 $\beta\:$Pictoris

The dusty and gaseous disk around the southern star $\beta\:$Pictoris ($\beta\:$Pic) (Smith & Terrile 1984) is considered today as a key example of a young planetary system in its early dynamical evolutionary phase, after the dissipation of the thick gaseous nebula (see reviews by Artymowicz 1997; Vidal-Madjar et al. 1998; Lagrange et al. 2000). It is known that the dusty particles viewed on the disk images are not a remnant of any primordial, more massive and opaque disk, but rather second generation material continuously replenished from inside the disk by larger, planetesimal-like bodies, either by slow evaporation (Lecavelier et al. 1996) or by collisions (Artymowicz 1997). Hence the presence of numerous planetesimal- or comet-like bodies within this disk is highly probable. The presence of larger bodies like planets is even suspected, either to explain the inner warp of the disk (Mouillet et al. 1997), or to explain the so-called Falling Evaporating Bodies (hereafter FEB) phenomenon that has been now monitored for years towards this star (see below).

The presence of kilometer-sized bodies was first introduced to explain the numerous spectral absorption events reported in the spectrum of $\beta\:$Pictoris. The survey of various spectral lines (Ca II, Mg II, Fe II, etc.) towards this star revealed that transient absorption features, usually redshifted, frequently appear in addition a deep central stable component. These additional features evolve within one day or even less (Boggess et al. 1991, Vidal-Madjar et al. 1994, Lagrange et al. 1996, and refs. therein).

These repeated spectral events have been successfully modeled as the signature of the evaporation of kilometer-sized bodies crossing the line of sight in the vicinity of the star, on star-grazing orbits (i.e. FEBs; Beust et al. 1990, 1996, 1998 [hereafter Paper XXV]).

A key issue concerning this scenario was the identification of a triggering dynamical mechanism capable to bring numerous bodies on star-grazing orbits, out of a Keplerian rotating disk on quasi-circular orbits. Various mechanisms were proposed, all of them involving the gravitational perturbations by at least one planet. As of yet, the most convincing mechanism was proposed by Beust & Morbidelli (1996, 2000): The star-grazers are generated from bodies initially orbiting the star on low eccentricity orbits, but trapped in a mean-motion resonance (the 4:1 one is by far the most powerful one) with a massive, Jovian-like planet orbiting the star on a moderately eccentric orbit ( $e'\sim 0.05$-0.1). Under such conditions, the trapped particles are able to become star-grazers within $\sim$10⁴ revolutions of the planet. This scenario reproduces fairly well from a statistical point of view the dynamical characteristics of the FEBs that were deduced from the observation of variable features.

1.2 Herbig stars

Attempts to find stars presenting variable spectral phenomena like $\beta\:$Pic were first carried out by Hobbs (1986) and Lagrange-Henri et al. (1990) for main sequence stars. As the FEB phenomenon is linked to the presence of circumstellar material, these searches focused on stars presenting an infrared excess or evidence for gaseous shells. A few interesting candidates showed up: HR10, HR2174 and HR6519 (51Oph). These stars were furthermore the subject of more intense investigations. In particular, HR10 (A main-sequence A2 type star) revealed peculiar spectral variability in its Ca II K line (Lagrange-Henri et al. 1990b), confirmed furthermore (Lecavelier et al. 1997; Welsh et al.1998). Whether these variable lines are due to FEBs events still needs to be confirmed.

Apart from these stars, redshifted (and less frequently blueshifted) features and spectral variations have been now observed for years towards several Herbig Ae/Be stars. From a strictly stellar point of view, Herbig Ae/Be stars are pre-main-sequence objects that may be regarded as $\beta\:$Pic precursors (see the reviews by Grady et al. 1996a, 2000). These observations are based on UV spectroscopic data gathered with IUE. In some cases, the Fe II lines as observed with IUE appeared asymmetrical with extended red wings (Grady et al. 1993, 1996b), indicating possible circumstellar components, sometimes redshifted by several hundreds of kms^-1. Over-ionized species indicating collisional ionization, like in the $\beta\:$Pic case are frequently observed. Moreover, whenever circumstellar components are observed towards doublet lines like Mg II k and h, the ratio between the two variable components often appears close to 1, indicating that the absorption is caused by opaque but clumpy material in front of the star (e.g. 51 Oph: Grady & Silvis 1993; Lecavalier et al. 1997). This characteristic is shared with $\beta\:$Pic, for which it was considered as a strong argument in favor of the FEB scenario. We shall now focus on three specific examples.

1.2.1 UX Orionis and similar stars

Spectacular spectral variations have been reported towards some Herbig stars in Na I and H$_\alpha$ lines. This is the case in particular for UX Ori (Grinin et al. 1994) and BF Ori (Grinin et al. 1996b; de Winter et al. 1999). In both cases, photometric variability is also reported, and interpreted as occultation by more or less opaque orbiting clouds; the spectroscopic variability itself was claimed to be due to a $\beta\:$Pic-like FEB phenomenon (Grinin et al. 1996a; de Winter et al. 1999). This is quite plausible, but the variability in Na I lines raises a crucial question: the lifetime of Na I atoms towards photoionization in the vicinity of a hot star like those under consideration does not exceed a few seconds. This shows first that Na I must be continuously produced, possibly by FEBs; but the lifetime of Na I is so short that the evaporated material may not have enough time to form a large enough cloud in front of the star, a necessary condition for generating deep absorption components. As a matter of fact, no FEB signature in Na I lines was ever detected towards $\beta\:$Pic. This question was investigated by Sorelli et al. (1996). Based on simulations, they show that Na I may be protected in a cloud of size comparable to that of the star, provided a hydrogen density cm^-3. This is valid for a solar chemical composition cloud. If one considers however metal-enriched gas that may result for example from the vaporization of a rocky body, the constraint may be lowered down to cm^-3. According to the authors, the FEB hypothesis is plausible. The required density remains nevertheless at least 100 times larger than the average density obtained in FEB simulations for $\beta\:$Pic (Beust et al. 1996). The amount of material in front of the line of sight for each of the UXOri observations corresponds according to Sorelli et al. (1996) to the complete vaporization of a 20 km sized body. If the process that brings the suspected FEBs towards the evaporating zone is gradual, their initial sizes should be considerably larger. The mere plausibility of the FEB hypothesis for these stars is thus still controversial. The least that can be stressed is that we are dealing here with a significantly more violent physics than towards $\beta\:$Pic.

1.2.2 AB Aurigae

AB Aur is a well-known and much studied, young (2-$4\,10^6\,$yr; van den Ancker 1998) Herbig Ae/Be star. This A0-Ve type star exhibits strong, variable wind components in various spectral lines. Its complex wind structure has been the subject of intense investigations in recent years (Catala & Kunasz 1987 [hereafter C87]; Bouret et al. 1997; Bouret & Catala 1998, 2000). It was modeled as a two-components spherically symmetric wind with a chromosphere, in addition to hot clumps. The derived, model-dependent, mass loss rate is $1.8\,10^{-8}~M_\odot\,$yr^-1.

ISO spectroscopy revealed the presence of circumstellar dust around this star (van den Ancker 2000; Bouwman et al. 2000) and an accretion disk detected by interferometric observations. HST/STIS imaging (Grady et al. 1999b) reveals the presence of circumstellar dust scattering the stellar light. As the scattered-light distribution is symmetric, Grady et al. (1999b) conclude that AB Aur is probably viewed close to pole-on, with an inclination with respect to the plane of the sky in any case less than $45^{\circ}$ (and more probably than 20- $25^{\circ}$). The fact that the star is viewed close to pole-on also explains why the wind can easily be modelled as spherically symmetric; this however allows us to question whether the wind is actually spherically symmetric or only axisymmetric, concentrated at high stellar latitudes.

Despite their complex structure which makes components identification difficult, the spectral lines of AB Aur exhibit short-term variations on their red wing that have been interpreted by Grady et al. (1999) as resulting from a FEB phenomenon.

1.2.3 HD100546

	Figure 1: Two IUE spectra of the k and h Mg II doublet ( $\lambda \lambda 2796.4$ and $2803.5\,$Å) in HD 100546, taken from Grady et al. 1997. These spectra were taken an March 7 and May 25, 1995. A clear variation is visible between the two epochs. The heliocentric velocity scales are relative to the two lines
Open with DEXTER

Perhaps the most convincing candidate for a $\beta\:$Pic-like phenomenon so far is HD 100546, a $\sim$10$^7\,$yr old B9 type Herbig star (van den Ancker 1998). It is an intermediate case between younger Herbig stars like AB Aur and a more evolved system like $\beta\:$Pic.

HD 100546 exhibits a very important infrared excess (IRAS), later on confirmed by ISO (Waelkens et al. 1996). Recently, a possible disk of circumstellar dust has been imaged around this star (Augereau et al. 2001; Pantin & Lagage 2001). Observations of the near-infrared spectrum of HD 100546 gathered with the ISO satellite (Malfait et al. 1998) reveal emission features characteristic for crystalline silicates, comparable to those of comets and also those observed towards $\beta\:$Pic itself (and also recently towards two younger systems; Sitko et al. 1999). Polarimetric variability has been reported and interpreted as resulting from light scattering by dust within the circumstellar environment (Clarke et al. 1999).

Striking spectral variable absorption components, redshifted by several hundreds of kms^-1 are reported in Grady et al. (1997) and Viera et al. (1999), towards several lines including those of over-ionized species, and also towards lines of more fragile elements like C I. We display here the Mg II IUE data taken from Grady et al. (1997) (Fig. 1), as these are data we will try to model. We note in Fig. 1 the presence in one of the two spectra of a deep variable feature appearing in both lines at $\sim$150 kms^-1 redshift with respect to the central absorption. 100-300 kms^-1 is indeed a typical velocity range for variable absorption features observed towards Herbig Ae/Be stars.

Note however that the characteristic time-scale of these variations is very poorly constrained by the observations presented in these papers. The spectra were indeed taken at different epochs separated by at least several months; one of the most striking characteristics of the FEB spectral variations towards $\beta\:$Pic is that they appear or disappear within one day or even less. The spectroscopic database available as of yet does not allow to test such variability time-scales.

Grady et al. (1997) show ratios equal to 1 between the variable components occurring in doublets, implying as for $\beta\:$Pic saturated clouds that do not mask the whole stellar disk. They derive filling factors of the stellar disk ranging between 40 and 50%. A rough estimate of the relative abundances of the metallic species in the variable lines shows consistency with a cometary-like chemical composition (with however more volatiles and magnesium than for $\beta\:$Pic), and the authors estimate the amount of material in front of the line of sight to the one corresponding to the entire evaporation of a 1 km sized body. Note however that the abundance estimate was performed assuming a close-box model which we showed not to be accurate in permanently refilled gaseous environments like FEB comas (Lagrange et al. 1998).

1.3 Modeling FEBs towards Herbig stars. Aim of this paper

The presence of a FEB-like phenomenon has then been proposed for many stars presenting spectral variations. However, such a global affirmation has to be taken with care. Obviously, a dedicated analysis of the detectability of spectral FEB signatures in the spectrum of a Herbig Ae/Be star is needed. As a matter of fact, their stellar environment is very different from that of $\beta\:$Pic. The main difference is that strong stellar winds originating from Herbig Ae/Be stars may have some influence on the ionic clouds escaped from the FEBs, and consequently affect their spectral signature. Modeling is needed, before claiming that a FEB scenario is plausible for accounting for the observations. This is the purpose of the present paper.

In the case of $\beta\:$Pic, we developed in the past years a code able to simulate the dynamics of the gaseous material escaped from an evaporating body in the vicinity of the star, together with the absorption components it generates when crossing the line of sight (Beust et al. 1990, 1996; Paper XXV). This theoretical tool helped us understanding the various aspects of the scenario, such as the ability to generate different kinds of variable features depending on the distance between the passing body and the star, or the behavior differences between various species.

We want now to adapt this modeling to the case of Herbig Ae/Be stars, in order to investigate in which conditions the FEB scenario can account for the observations. This parametric study is intended to generally apply for any Herbig star. We shall apply afterwards this study to AB Aur and HD 100546.

In Sect. 2, we describe into more details the model we assume for Herbig stars. We point out the differences we need to assume between such a star and $\beta\:$Pic, in particular the role of the stellar wind. In Sect. 3, we describe the simulations relevant to FEBs around Herbig stars. In Sect. 4, we discuss our results and their application to the quoted examples of Herbig stars. Our conclusions are presented in Sect. 5.

2 A model for a Herbig star. Comparison with $\beta\:$Pic

2.1 Basic features

The input stellar parameters needed for our simulation code are the following: the mass M and the radius R of the star, its effective temperature $T_{{\rm eff}}$, and a calibrated spectrum ranging from the UV to the infrared. The mass of the star is required to compute the dynamics of particles. The calibrated spectrum is the key stellar parameter. It is needed to compute the evaporation of dust particles, their dynamics, the dynamics of metallic ions, and the synthetic spectrum when some ions cross the line of sight.

	Figure 2: The radiation pressure ($\beta$ value) acting on dust grains as a function of the grain radius, for $\beta\:$Pic (thin lines), a B9 type star like HD 100546 (thick lines), and the Sun (dashed lines). For each star, two curves are given, one corresponding to grains with zero porosity, and another one corresponding to grains with 0.95 porosity
Open with DEXTER

Figure 3: Sublimation region for dusty "astronomical silicates'' grains as a function of their size, in the vicinity of $\beta\:$Pic, of a B9 type star and of the Sun. In each case, the sublimation zone is basically delimited by the drawn grey band: the inner edge of this region corresponds to the distance where the lifetime of the grains towards sublimation is 0.1day, while the outer edge corresponds to a lifetime of 10 days. The dashed line indicates a lifetime of 1 day

Open with DEXTER

Along the path of the evaporating body, the simulation releases continuously fresh particles close to the nucleus. These new particles are first assumed to be dust particles chosen with a convenient Hanner distribution function (see Paper XXV for details). At each time step, the stellar spectrum is used to compute the equilibrium temperature of the dust grains via Mie theory (depending on their stellar distance), and the strength of the stellar radiation acting on each of them. The radiation pressure itself is described as usual by a its ratio $\beta$ to stellar gravity; $\beta$ does not depend on the distance to the star, but for dust grains it depends on their size. Also the evaporation rate (function of the equilibrium temperature) is computed for each grain, together with its radius decrease within the considered time-step. The $\beta$ ratio is used to compute the motion of the grain during the same time-step.

We model our dust particles as "astronomical silicates'', assuming the optical constants of Draine & Lee (1984), revised in Laor & Draine (1993). Once a given dust particle has reached zero radius, it is then assumed to represent metallic ions. Its dynamics is then computed taking into account i) the radiation pressure (i.e. the $\beta$ ratio) acting on it, eventually reduced by a shielding effect from the other ions located closer to the star, and ii) the drag force by the neutral gas escaped from the nucleus. This is described extensively in Beust et al. (1990, 1996). Finally, the spectral absorption of the ions passing in front of the star is computed, and a synthetic spectrum is produced.

As a template model for a Herbig Ae/Be, we assume a B9 type star with $M=2.4~M_\odot$, $R=1.75\,R_{\odot}$, $T_{{\rm eff}}=10\,500\,$K, $\log g=4.0$, $v\sin i=250\,$kms^-1and compute an ATLAS9 Kurucz synthetic spectrum (Kurucz 1979) relevant for these parameters. These parameters correspond actually to a model for HD 100546 (van den Ancker 1998), but may be considered as typical for any Herbig Ae/Be star. For AB Aur (A0 type star), a more exact set of parameters is (van den Ancker 1998) $R=2.5\,R_{\odot}$, $T_{{\rm eff}}=10\,000\,$K, $\log g=4.1$, $M=2.4~M_\odot$. Assuming the latter set of parameters instead of the former appears to only slightly affect the numerical results of this study and does not change the conclusions we derive below. Consequently, we will not present separate results for both sets of parameters. We will only present those corresponding to the B9 type star model, considering they are typical for Herbig Ae/Be star in general.

Figure 2 plots the $\beta$ value as a function of grain size, for two values of grain porosity: 0 (a solid grain), and 0.95 (a more realistic cometary grain value), and for three stellar cases: a B9 star, $\beta\:$Pic, and the Sun. In the case of $\beta\:$Pic, the parameters assumed are (Crifo et al. 1997): $M=1.7~M_\odot$, $T_{{\rm eff}}=8200\,$K, $R=1.5\,R_\odot$. The spectrum is a Kurucz model computed with the quoted value for $T_{{\rm eff}}$ and $\log g=4.25$.

The curve corresponding to $\beta\:$Pic and porosity 0 is to be compared to that published by Artymowicz (1988). We first see that the shapes of all curves are similar: smaller grains undergo a much stronger radiation pressure than larger ones, and tend to be blown away by the star ($\beta>1$). For smaller grains, the radiation pressure tends to be $\sim$5 times larger in a B9 star environment than around $\beta\:$Pic, while for larger grains, the difference tends to decrease.

In all cases, grains with porosity 0.95 undergo a radiation pressure at least one order of magnitude larger than the grains with no porosity. We will assume in the numerical simulations that the porosity of the grains is 0.95 as it is a a more realistic cometary value.

A more striking difference shows up in Fig. 3. Here the region of dust sublimation in the various stellar environments is shown as a function of the grain radius. More precisely, the grey zone corresponds to the stellar distance range where the lifetime of the grains lies between 0.1day and 10days. Most of the sublimation of the dust escaped from FEBs is expected to occur in this area.

We first see, as noted in Beust et al. (1998), that smaller grains tend to evaporate further away from the star(s) than larger ones. In the case of a B9 star, the grains evaporate at larger distances that around $\beta\:$Pic, and that around the Sun, the refractory grains are expected to evaporate only in the immediate vicinity of the star. This is a direct consequence of the difference in effective temperature between the various stars. The difference is more important for smaller grains (radius $1\,\mu$m). Note that the FEBs are supposed to release grains of all sizes, with a given size distribution. They can generate transient absorption events as soon as some of these grains begin to evaporate. Hence FEBs are expected to be observable around $\beta\:$Pic as soon as they get a periastron 0.5AU, while in the environment of a B9 star, this limit jumps up to $\sim$2AU, i.e., 4 times further away. This is a major difference that tends to indicate that active FEBs around Herbig stars can be more frequent than around $\beta\:$Pic, but whether these FEBs can generate detectable signatures still needs to be investigated. Note conversely that the dust sublimation region around the Sun does not exceed 0.1AU.

 

 

Table 1: Compared $\beta$ values for various species in the environment of $\beta\:$Pic and that of a B9 type star. The values for $\beta\:$Pic correspond to those listed in Lagrange et al. (1998)

Element	$\beta$	$\beta$	Element	$\beta$	$\beta$
	($\beta\:$Pic)	(B9 star)		($\beta\:$Pic)	(B9 star)
Ca II	35.1	318.9	Fe I	19.1	90.9
Mg II	5.0	84.4	Fe II	4.9	46.4
Mn II	11.0	62.2	Ni II	1.39	15.0
Al II	0.56	132.9	Zn II	6.1	47.64
Al III	15.6	108.1	Cr II	3.4	23.9
Si I	10.1	136.1	S I	1.52	63.7
Si II	0.15	13.9	C I	0.095	41.4
H I	$\sim$0.001	0.43	Na I	42.8	477.1

Let us now compare in both stellar environments the behavior of the metallic ions that are produced by the dust evaporation. Basically the main parameter is the radiation pressure, which is computed for a given element by the adding the contribution of all its spectral resonance lines (see Lagrange et al. 1998, 1996 for details), and described as for dust grains with a $\beta$ ratio. $\beta$values for $\beta\:$Pic a B9 star are given for some important species in Table 1. The B9 star values are significantly larger than the corresponding $\beta\:$Pic ones. This is clearly due to the more important stellar flux of a B9 star. The difference is even more striking for some elements like C I or H I. In that case, most of the contributing spectral lines are located fairly far in the UV spectrum ( $\lambda<1700$Å). Downwards these wavelengths, the $\beta\:$Pic flux drops sharply, while this is not the case for a hotter star.

In any case, we can stress that two identical FEBs, with similar evaporation rates, located at the same stellar distance, one in the vicinity of $\beta\:$Pic and the other one in the vicinity of a Herbig star, will probably generate ionic clouds with different sizes. The size of a ionic cloud in front of the star is roughly fixed by the distance to the nucleus where the radiation pressure exactly balances the drag force from the surrounding neutral medium (see Beust et al. 1990; Papier XXV). In the B9 star case, as the radiation pressure is significantly larger than for $\beta\:$Pic, we expect this balance to occur closer to the nucleus, and consequently the resulting ionic clouds to be smaller. As the projected size of the cloud directly controls the maximum depth of the variable spectral component it can generate, we thus expect the spectral components generated by identical FEBs to be fainter in the case of a Herbig star than for $\beta\:$Pic, and therefore more difficult to observe.

2.2 Evaporation as a function of stellar distance

Just like the refractory material is able to evaporate further away from the star for earlier spectral types, we expect comet-like bodies to start releasing volatile material, and consequently dust grains, at larger stellar distance in the same case. This will basically depend on the sublimation of water ice, the major volatile component. For Solar System comets, the evaporation rate indeed depends on the stellar distance, but also on the nature of the body, on its rotation, etc. Seasonal effects are also observed. Assuming a simple isothermal nucleus model, it is nevertheless possible to derive a simple description that depends only on the heliocentric distance, and that provides a convenient first-order fit of the observations. Rickman (1991) gives the following law for the H₂O sublimation flux as a function of heliocentric distance:

	Figure 4: Surface sublimation rate for a large ice sphere as a function of stellar distance, in the environment of i) $\beta\:$Pic (thin solid line), ii) a B9 type star (thick solid line), and iii) the Sun (dashed line). In addition, we have plotted (dotted line) the empirical function g(r), once rescaled to the Solar curve
Open with DEXTER

 

$$0.111262\left(\frac{r}{2.808\,\mbox{AU}}\right)^{-2.15} \;$$ g(r) =

$$\qquad\mbox{}\times\left[ 1+\left(\frac{r}{2.808\,\mbox{AU}}\right)^{5.093}\right]^{-4.6142}.$$ (1)

This dimensionless function is normalized to g=1 at $r=1\,$AU. For stars hotter than the Sun, we expect the behavior to be somewhat different. This may be roughly evaluated, computing the evaporation rate of ice from basic Mie theory in each stellar environment. We already described this formalism in Beust et al. (1998). The mass sublimation rate per unit of surface of the body is given by (Lamy 1974):

$$\frac{{\rm d}E}{{\rm d}t}=0.7p(T)\sqrt{\frac{m}{2\pi kT}},$$ (2)

where T is the equilibrium temperature of the body, p(T) is the vapor pressure, k is Boltzmann's constant, m is the molecular weight of the ice.

The equilibrium temperature itself is derived from the energy balance between the stellar radiation flux, the thermal energy re-radiated by the body, and the sublimation energy (see Beust et al. 1998). The absorption efficiency is obtained via Mie theory, but for a large body as we assume here (we use a typical size of 1 km), we are well inside the geometrical optics regime, so that the equilibrium temperature, and consequently ${\rm d}E/{\rm d}t$ (as long as it may be assumed homogeneous) are virtually independent from the size of the body. The latent heat for ice sublimation and the vapor pressure are taken from Lamy (1974), and are listed in Paper XXV.

Figure 4 gives the computed value of ${\rm d}E/{\rm d}t$ as a function of stellar distance, for the three stellar environments of $\beta\:$Pic, our template B9 type star, and the Sun. In addition, we have superimposed (dotted line) the g(r)profile, once rescaled to our computed Solar curve. The agreement is very good.

We first note that, as expected for a given stellar distance, the hotter the star, the larger the sublimation rate. Basically for stellar distances less than 1AU, ${\rm d}E/{\rm d}t$ is 11 times larger in the $\beta\:$Pic environment that around the Sun, and it is 3 times larger around a B9 star than around $\beta\:$Pic. We also note that the sublimation of ice begins at much larger stellar distances for hotter stars. For Solar comets, the sublimation rate drops sharply further than $\sim$2AU from the Sun. Around $\beta\:$Pic, this limits jumps up to $\sim$10AU, and for a B9 star, it reaches $\sim$15AU. This means that we should expect FEBs around Herbig stars i) to have more important outgassing rates than $\beta\:$Pic one, and ii) to start activity further away from the star.

The actual sublimation rate of a given body may be somewhat different than predicted in Fig. 4. It can depend on the surface characteristics of the body, such as for instance its Bond albedo. Moreover, this is only a rate per unit surface. The whole production rate of a given FEB should in fact also depend on its size. We may nevertheless stress that two identical bodies located at the same stellar distance from $\beta\:$Pic and a B9 star respectively, should present outgassing rates within a factor $\sim$3, provided they are close enough to their stars ( 1AU).

In our simulation runs, our strategy concerning production rates will thus be the following: once a FEB orbit is chosen, we give as input the dust and gas production rates, taking care that the dust/gas mass production rate is of order unity. We assume these values to apply at periastron, and we scale the production rates at larger stellar distances according to the ${\rm d}E/{\rm d}t$ law plotted in Fig. 4, relevant for the stellar environment under consideration. The $\beta\:$Pic spectral events are well simulated assuming typical dust and gas production rates of $\sim$ $3.5\,10^{33}$mol.s^-1 and $\sim$$5\,10^7\,$kms^-1respectively, for a body passing at 0.2AU. These values may be considered as accurate within a factor 3 or 4 only. We will nevertheless assume them as a starting point of our analysis in the environment of a B9type star, first multiplying them by 3 according to Fig. 4.

2.3 Stellar wind models

Perhaps the most striking difference between the environment of a main-sequence A-type star like $\beta\:$Pic and a Herbig Ae/Be star is the presence of a strong stellar wind. No stellar wind has been detected so far originating from $\beta\:$Pic. From a theoretical point of view, A-type stars are expected to have radiatively driven winds made only of metals (those which indeed suffer most of the radiation pressure), with a mass-loss rate less than $10^{-16}~M_\odot$yr^-1 (Babel 1995). Conversely, towards pre-main-sequence objects like T Tauri and Herbig Ae/Be stars, there are spectroscopic evidences for strong stellar winds characterized by terminal velocities ranging between 100 and 500 kms^-1, with global mass-loss rates of $\dot{M}\sim 10^{-9}$- $10^{-7}~M_\odot$yr^-1 (see reviews by Bertout 1989; Catala 1989).

Focusing now to our two example stars, we note first that the stellar wind of AB Aur has been well studied and modeled. Bouret & Catala (1998) derive a mass loss rate of $(1.8\pm0.3)\,10^{-8}~M_\odot$yr^-1, which falls in the quoted range. There is so far no specific measurement on the mass loss rate of HD 100546, but Viera et al. (1999), on the basis of the analysis of the H$_\alpha$ profile of the star, conclude to the presence of a stable wind at velocity $\sim$400 kms^-1, in addition to more or less variable accretion episodes.

Obviously the presence of a strong stellar wind should drastically affect the coma of the suspected FEBs, and modify their spectroscopic signatures with respect to what should be expected without any noticeable wind. Hence any proper modeling of the FEB scenario around Herbig Ae/Be stars must take the interaction with the stellar wind into account.

It is first necessary to state which model we should assume for a stellar wind of such a star. This is in fact not straightforward as the origin of Herbig Ae/Be winds is still unclear, if not controversial. Basically, three major mechanisms are able to push a stellar wind:

Radiation pressure

can blow away gaseous material from the stellar atmosphere and generate a wind. In fact, this is only relevant for O-type stars (or earlier) (Castor et al. 1975), as only those stars are able to push away hydrogen;

The thermal expansion

of a high-temperature gas is the classical mechanism that generates winds for Sun-like stars. This process is nevertheless efficient only if there is a hot corona (T>10⁶K) within a few stellar radii from the stellar surface. The classical model for this kind of wind is described by Parker (1958). The high coronal temperature is usually thought to be linked to the presence of an extended convection zone under the stellar surface. Hence this model is probably not relevant for main-sequence stars of early types such as $\beta\:$Pic, because of their lack of significant convection zone under the surface. Herbig Ae/Be stars are however not main-sequence stars, and may thus have large convection zones. Indeed, X-ray emission has been observed towards many such stars (Damiani et al. 1994).

Hydromagnetic waves

can also generate stellar winds. If the star possesses a large-scale fossil magnetic field, the turbulent motion at the base of the field lines will induce oscillations that will propagate as Alfvén waves. If these waves are not damped, they may accelerate ionized particles located around the star and generate a strong wind (see a description by Strafella et al. 1998). This mechanism can be very important for Herbig Ae/Be stars, but it is still unclear. The description of such a wind is nevertheless far from being straightforward. The geometry of the flow may even not assume spherical symmetry. In our parametric study, we shall limit ourselves here for simplicity to spherically symmetrical winds. Basically, we will assume two different kinds of models.

The first one is a solar-like hot isothermal ($T=10^6\,$K) coronal wind (Parker 1958). Once the temperature is fixed, this model depends only on one single parameter, namely the mass loss rate $\dot{M}$. The wind velocity profile has to pass through the sonic critical point, which location is fixed by the value of $\dot{M}$. The velocity law v(r) then assumes the form

$$v(r)\exp\left[-\left(\frac{v(r)}{c_{{\rm s}}}\right)^2\right] ={\rm Cst.}\;,$$ (3)

where $c_{{\rm s}}$ is the sound speed.

The second model will assume is the semi-empirical model introduced by C87 for Herbig Ae/Be stellar winds. They assume the following temperature profile T(r), depending on the four parameters $(T_0,T_{{\rm max}}, \Delta_1,\Delta_2)$:

$$T(r)=\left\{\begin{array}{rcl} \multicolumn{3}{l}{\dy T_0+(T_... ...3.5truecm} & \mbox{if} & r\geq R_{{\rm ch}}\end{array}\right.,$$ (4)

with

$$R_{{\rm ch}}=R_*+\frac{\Delta_1}{2\sqrt{\ln 2}}\sqrt{\ln\left( \frac{T_{{\rm max}}-T_0}{T_{{\rm eff}}-T_0}\right)}\;.$$ (5)

This temperature law is well suited to mimic a chromospheric temperature bump above the photosphere. The maximum temperature $T_{{\rm max}}$is reached at $r=R_{{\rm ch}}$, which value is fixed as a function of the other parameters via Eq. (5). This just ensures that the temperature at the level of the stellar photosphere (r=R_*) equals the effective temperature of the star ( $T_{{\rm eff}}$). T₀ is the asymptotic temperature value at infinity. $\Delta_1$and $\Delta_2$ control respectively the position of the temperature maximum ( $R_{{\rm ch}}$) and the extension of the chromospheric region. In practice, the asymptotic temperature T₀ is almost reached at $r\simeq 2\Delta_2$, so that after the chromosphere, the wind is virtually isothermal.

Using this model to fit the chromospheric wind lines they observed, Bouret & Catala (1998) gave estimates for all these parameters for four typical Herbig Ae/Be stars, including AB Aur. Schematically, the typical values they get are , , , and , depending on the peculiar star under consideration. In the following we will assume in our tests the values quoted by Bouret & Catala (1998) for AB Aur, namely $T_{{\rm max}}=17\,000\,{\rm K}$, $T_0=4800\,{\rm K}$, $\Delta_1/R_*=0.17$ and $\Delta_2/R_*=1.5$. Note that in any case, the temperature remains considerably less than the one we assume in the coronal wind model ($10^6\,$K).

	Figure 5: The two velocity profiles assumed in our wind models. The dashed line is the coronal wind profile, and the solid line is a piecewise linear function taken from C87. The velocity gradient of the latter is much steeper in the first stellar radii above the surface, but the terminal velocity is less
Open with DEXTER

	Figure 6: Schematic representation of the global interaction between the stellar wind and the cometary atmosphere
Open with DEXTER

If the physical process that causes the wind expansion was clearly identified, then it would be possible, as for the coronal wind model, to derive the velocity profile of the wind from the usual hydrodynamical conservation equations. As this is not the case, C87 assume empirical piecewise linear velocity laws, each of them starting with negligible velocity at stellar surface, a steep velocity gradient within the chromosphere, and a terminal velocity of 300 kms^-1. We chose to assume the velocity law labeled "V₁'' in C87, as it represents the best fit for AB Aur. This profile and the coronal wind one are displayed in Fig. 5.

Once the velocity profile is fixed, the density profile $\rho(r)$is obtained from the continuity equation

$$\dot{M}=4\pi r^2\rho(r) v(r)\;.$$ (6)

2.4 Interaction between the wind and the FEBs

In presence of a strong stellar wind, we expect the dynamics of the gaseous material escaped from the nucleus of the FEBs to be drastically affected. Similarly, the flow of the wind around the FEB nucleus should also be affected, like for the solar wind around comets.

The interaction of comets with the solar wind has been described by several authors (see e.g. Galeev et al. 1985; Flammer 1991; Flammer et al. 1997); it is sketched in Fig. 6. Basically, the solar wind is affected because the gaseous material escaped from the nucleus tends to be ionized by the ultraviolet radiation field or by charge exchange interaction with the wind particles. The new ions produced are then incorporated to the wind. This mass loading process decelerates the wind in front of the comet. If the deceleration is strong enough (which is usually the case), then a bow shock forms. The supersonic wind becomes subsonic after the shock. Closer to the nucleus, the wind is then stopped along a stagnation surface (or tangential discontinuity). Inside this surface we find only cometary plasma, which undergoes itself an inner shock before reaching the tangential discontinuity. This picture was confirmed by in-situ measurements. For typical solar comets at $\sim$1AU from the Sun, the first shock forms at $\sim$5 10⁵ km from the nucleus, while the wind stagnation takes place at $\sim$10⁴ km.

In our description, we will assume that the flows is axisymmetric around the line comet - star. Following Krankowsky (1991), we will assume for simplicity that the gaseous material escaped from the nucleus is made of water vapor with a number abundance of 80%, and that the 20% left consist of carbon monoxide (CO). In the environment of a Herbig Ae/Be star, we expect the water molecules to be quickly photo-dissociated. The resulting H and O atoms should not be subject to stellar photoionization (the star is not hot enough), so that they do not take part to the interaction with the wind. Conversely, CO molecules are able to resist a long time enough to allow the formation of CO⁺ ions which may interact with the stellar wind. As a matter of fact, water photodissociation and CO photoionization were listed by Huebner et al. (1991) as the most important chemical reactions expected in the coma of a comet.

We also assume that inside the tangential discontinuity there is no wind, and that the neutral flow is not perturbed. The details of the dynamics of the wind flow before and after the bow shock are described in the appendix, as well as the shape of the shock and of the tangential discontinuity. The wind is assumed to come from infinity in a plan-parallel geometry; we neglect the transversal motion of the wind around the nucleus. The basic parameters are the velocity and the density of the wind at infinity, which depend on the distance between the comet and the star, via the wind model assumed; the volatile production rate Q, outflow velocity $v_{{\rm e}}$, ionization time $\tau$, and recombination rate $\alpha$. For solar comets, $v_{{\rm e}}$ is usually about 1 kms^-1 and $\tau\simeq 10^6\,\mbox{s}\times(d/\mbox{1\,AU})^2$(Flammer et al. 1997), where d is the distance to the star. We will assume the same values. The recombination rate $\alpha$ does not depend sharply on the temperature. Following Gombosi et al. (1983), we assumed $\alpha=10^{-6}\,\mbox{cm}^3\,\mbox{s}^{-1}$. From these input parameters, we are able to derive the characteristics of the wind flow at any point within the cometary coma.

The present description holds in fact for a comet which does not move very fast with respect to the wind. This is justified for Solar comets at 1AU from the Sun which have an orbital velocity $\approx$40 kms^-1 while the Solar wind velocity is more than 400 kms^-1. This is far less justified for FEBs passing at periastron at a few stellar radii from the stellar surface: closer to the star, the orbital velocity is larger while the wind velocity is less. At $\sim$10 or 20 stellar radii from the surface (depending on the wind model and of the stellar parameters), both velocities turn out to be comparable.

A convenient way to (partially) solve this problem is to take into account the velocity of the comet when computing the velocity of the wind at infinity $\vec{u}_{\infty}$ (see Appendix). If the comet is located at distance $r_{{\rm c}}$ from the star, where the wind velocity is $\vec{u}_{{\rm wind}}$, and moves with velocity $\vec{u}_{{\rm comet}}$, we may just assume (vectorially)

$$\vec{u}_{\infty}= \vec{u}_{{\rm wind}}-\vec{u}_{{\rm comet}}\;.$$ (7)

Once this is assumed the same formalism applies, the main difference being that, as the apparent incoming velocity of the wind $\vec{u}_{\infty}$ is no longer necessarily aligned with the direction comet - axis, the bow shock and the tangential discontinuity now appear more or less inclined with respect to that axis, depending on the relative velocities.

We want to point out here that even with this refinement, the present description is probably a fairly crude approximation of the stellar wind flow around a more or less star-grazing comet. If the comet is close enough to the star, the approximation of a plan-parallel wind flow coming from infinity is very bad. One should consider the actual shape of the (spherical) wind expansion around the star, and compute how it is modified in presence of a FEB. Moreover, as the FEB passes near perihelion, it turns quickly around the star. This causes centrifugal terms to appear when the dynamics of the flow is computed in a referential frame bound to the nucleus. These effects are indeed taken into account in our simulation code when we compute the dynamics of individual metallic ions escaped from the FEB (Beust et al. 1996). Taking them nevertheless into a 3-dimensional hydrodynamical computation of the wind flow around the star-grazing comet is far beyond the scope of the present work, as this would require a specific study, probably highly computing time consuming. The purpose of this paper is only to seek the first-order effects of the introduction of a strong stellar wind. Moreover, our knowledge of the structure and the geometry of Herbig Ae/Be winds is so rough that getting into more details in the description of the interaction with the cometary medium would probably not be relevant. We decided however arbitrarily to consider that if the comet is located at distance d from the star, the zone of perturbed wind, as computed from the present model should not extend further away than 2d from the star. Beyond this limit, we assume that the wind recovers its unperturbed geometry.

Now, we need to describe the interaction between a wind flow (perturbed or not) and the individual metallic ions escaped from the FEB, and supposed to be responsible for the transient spectral events observed. We recall that once escaped from the FEB nucleus, the metallic ions are basically subject to two major conflicting actions, i) the radiation pressure that tends to blow them away in the anti-star direction, and ii) the drag (or collision) force due to the surrounding volatile medium that tends to retain them around the nucleus. Now we must add an interaction force with the wind, i.e. a ionized flow. We had already described this interaction in Beust et al. (1989). It results from classical Coulomb scattering theory. Let us consider one ion of mass m and charge q, moving at velocity v in a field of protons (charge e). After one collision, the impulsion change in the direction parallel to the unperturbed velocity reads

$$\delta p_{\parallel}=\mu v(\cos\chi-1)$$ (8)

where $\mu$ is the reduced mass, and $\chi$ is the deflection angle given by:

$$\tan\frac{\chi}{2}=\frac{1}{4\pi\epsilon_{0}}\frac{q{\rm e}}{\mu}\frac{1}{bv^{2}} \equiv\frac{C}{bv^{2}},$$ (9)

with b denoting the impact parameter. Averaging over all collision occurrences, we derive the drag force

$$f=nv\int_{0}^{\lambda_{{\rm D}}}2\pi b\delta p_{\parallel}\,{\rm d}b,$$ (10)

where n is the number density of the ionized medium, and where the integral has been classically cut towards large impact parameters to the Debye screening length $\lambda_{{\rm D}}$ of the plasma. After integration, we get

$$f=-2\pi n\mu\frac{C^2}{v^2} \ln\left(\frac{\lambda_{{\rm D}}^{2}v^{4}}{C^{2}}+1\right)$$ (11)

where C is defined in Eq. (9).

   
3 Simulations

Figure 7: Schematic plots of two simulation outputs, for a FEB passing at $\sim$0.1AU from the star ($q=0.04\,$AU) and a wind mass loss rate $\dot{M}=10^{-8}~M_\odot$yr-1, assuming a) a Parker (1958) wind structure, and b) a C87 one. In both cases, the upper plot shows a 2-D projection (the simulation is actually 3-D) in the orbital plane (which is assumed to contain the line of sight) of the Mg II cloud that expands around and behind the evaporating FEB nucleus, together with a small portion of the orbit and the star itself. The line of sight is sketched as two parallel dashed lines. The lower plots show the synthetic Mg II h and k bottom line profiles corresponding to this situation. In each case, the dotted profile superimposed to the actual one is the unperturbed profile that is to be observed when no FEB is present, so that the difference between the two profiles is due to the FEB. Note that contrary to similar plots concerning $\beta\:$Pic that were presented in earlier papers, we have not added here any stable circumstellar absorption centered at the star velocity, as this is more specific to $\beta\:$Pic itself. Moreover we have arbitrarily set the radial velocity of the star to 0, so that the heliocentric velocities coincide with velocities relative to the star. Due to the stellar wind action, the Mg II clouds appear in fact as a very thin tail which is unable to cover a significant portion of the stellar surface, so that the resulting spectral signatures are virtually undetectable. For clarity, we did not display here on the plot the location of the shocks, because they are located very close to the nucleus. Conversely, they are shown in the enlargement in Fig. 8

Open with DEXTER

3.1 The choice of the parameters

The present formalism has been included into our simulation code for FEBs. We choose to present here simulation results for Mg II ions, as these may be easily compared to available data for HD 100546 and AB Aur.

In the FEB scenario, the bulk redshift velocity v of a variable absorption corresponds to the projection onto the line of sight of the orbital velocity of the corresponding FEB. This results from the fact that most of the absorption is due to the ionic coma that surrounds the nucleus, and not to the tail. Assuming a Keplerian formalism, we easily derive this velocity v (Beust et al. 1996) as a function of the periastron value of the FEB (or its stellar distance d when crossing the line of sight), its eccentricity e and its longitude of periastron $\varpi$with respect to the line of sight:

$$\sqrt{\frac{GM_*}{d}}\frac{{\rm e}\sin\varpi}{\sqrt{1+{\rm e}\cos\varpi}} \;=\;\sqrt{\frac{GM_*}{q(1+{\rm e})}}{\rm e}\sin\varpi$$ v =

$$\simeq \frac{46\,\mbox{km\,s}^{-1}}{\sqrt{d/\mbox{1\,AU}}} \frac{{\rm e}... ...\,s}^{-1}} {\sqrt{q/\mbox{1\,AU}}}\frac{{\rm e}\sin\varpi}{\sqrt{1+{\rm e}}}\;,$$ (12)

where we have assumed $M_*=2.4~M_\odot$; e is supposed to be close to 1 as the orbit is star-grazing. Requiring $v\simeq 150\,$kms^-1 as seen from Fig. 1 turns out to be very constraining for these parameters. For any given periastron value q, the maximum velocity is reached for $\varpi=90^\circ$, and even in that case, $v\geq150\,$kms^-1can only be reached with AU, which is very close to the star. As a matter of fact, we had already shown that in the case of $\beta\:$Pic, the highest velocity events sometimes observed (High Velocity Features or HVFs) at kms^-1 could only be due to FEBs passing very close to the star with high values of $\varpi$ (Beust et al. 1998). Conversely the more frequent Low Velocity Features (LVFs) appearing at $\sim$20 kms^-1 in the $\beta\:$Pic spectrum are well reproduced with FEBs passing somewhat further from the star ($\sim$0.2AU) more closely before periastron ( $\varpi\simeq 20^\circ$).

In order to reproduce $v\simeq 150\,$kms^-1 in the B9 star environment, we must then assume a FEB passing closer to the star than 0.05AU and with orbital orientation more or less perpendicular ( $\varpi\simeq90^\circ$) to the line of sight. If we assume an even smaller periastron value, then the constraint on $\varpi$ is weaker, but this does not change anything to the results presented below. We will then take in all the following simulations a typical FEB orbit with $q=0.04\,$AU and $\varpi=90^\circ$.

As mentioned above, in any simulation run, we must fix the outgassing and dust production rates at periastron. We start from the values fitted at 0.2AU in the $\beta\:$Pic environment, first multiplying them by 3 according to Fig. 4 to take into account the difference between the stellar environments, and multiplying them once again by (0.2/0.04)²=25 to take into account the closer periastron value. This yields a volatile outgassing rate of $2.6\,10^{35}$mol.s^-1 and a dust production rate of $3.7\,10^{9}~$kgs^-1. These will constitute our reference values, intended to be valid for a FEB identical to those assumed around $\beta\:$Pic, but located at 0.04AU from a B9 type star. Remember that all throughout the simulation run, the instantaneous production rates of the body under consideration are computed according to Fig. 4 with respect to the values at periastron given as input.

	Figure 8: A close up view of the cometary head zone corresponding to Fig. 7b. The location of the bow shock and the tangential discontinuity are now displayed on this enlargement. The flow of Mg II ions is severely constrained by the interaction with the wind
Open with DEXTER

3.2 Simulation results

Figures 7a and b show the result of simulations for a FEB for which we have assumed these orbital and outgassing parameters. In Fig. 7a the Parker (1958) wind model was assumed, while in Fig. 7b, it is the C87 one. In both cases, the mass loss rate was fixed to $10^{-8}~M_\odot$yr^-1, as a typical value for a Herbig Ae/Be star.

Figure 9: A simulation of the same FEB than in Fig. 7b, but in the $\beta\:$Pic environment, i.e., with a stellar wind mass loss rate $\dot{M}=10^{-16}M_\odot$yr-1, and production rates that have been lowered by a factor 11, according to Fig. 4. The wind is so weak that the shock surfaces appear to be pushed away up to the star itself. In this case of course, our treatment of the interaction between the wind and the FEB coma is no longer appropriate, but from a qualitative point of view, it may be stressed that here, contrary to the Herbig Ae/Be case, the dynamics of the Mg II ions is fully controlled by the radiation pressure and not by the wind. Consequently the Mg II cloud is large and a big variable features appears in the spectrum

Open with DEXTER

	Figure 10: Same as Fig. 7b, but now with a stellar wind mass loss rate $\dot{M}=10^{-11}M_\odot$yr-1. Contrary to Fig. 7, the bow shock and the tangential discontinuity have been here sketched here, as they appear now further away from the nucleus. Detectable transient absorption components due to the FEB are now observable
Open with DEXTER

	Figure 11: The spectral signature in the Mg II k line of the same FEB as in Figs. 7 and 10, for different stellar wind mass loss values
Open with DEXTER

We first note that there are very few differences between the two plots displayed, though the wind models are different. This is due to the fact that the structure of the wind flow around the nucleus (i.e., the location of the shocks) is essentially a function of the impulsion flux of the wind $\rho(r)v(r)$ only, which is itself directly related to the mass loss rate $\dot{M}$ by the continuity Eq. (6). Thus, the global action of the wind on the FEB at a given stellar distance is to first order controlled by the mass loss rate, irrespective of the wind model assumed. Note that this is only valid for a spherically symmetric wind. In the following, we will then not distinguish furthermore between the results for the Parker (1958) wind model and those for the C87 one, as in any case we checked that both were very similar. We will only present results corresponding to the C87 model, as this represents a more realistic wind model for a Herbig Ae/Be star.

Second, we note that the ionic cloud around the FEB is very thin. This is in fact due to interaction with the wind, and shows up clearly in Fig. 8, which is an enlargement of the cometary head zone for the Fig. 7b case (the C87 wind case), and where the location of the shock surfaces with the wind (bow shock and tangential discontinuity) has been superimposed to the plot for clarity. Note the inclination of these surfaces with respect to the star-FEB axis. This is due to the velocity of the FEB, which is comparable to the wind velocity. Now, the dynamics of the Mg II ions appears obviously strongly affected by the interaction with the wind. They indeed hardly expand further away than the tangential discontinuity (recall that there is no wind inside this area) before flowing along it. In fact, the size of the cloud appears in the present case fully controlled by the wind.

The major consequence of the small size of the ionic cloud is that the transient absorption event generated by the FEB is not detectable. Independently from line saturation, the cloud cannot absorb more flux than its relative projected size onto the stellar surface. Here this size is so small that the variable absorption is unobservable. This is the striking difference with the situation for identical FEBs orbiting $\beta\:$Pic. As a matter of fact, if we simulate the same FEB around $\beta\:$Pic, and adding a small wind of $10^{-16}~M_\odot$yr^-1, the wind appears to have virtually no role, as the shocks appear much further away from the nucleus than the typical size of the ionic cloud as it is sculpted by radiation pressure. In the $\beta\:$Pic case, a strong transient event can be observed (Fig. 9).

The net effect of a strong stellar wind on the FEB phenomenon is thus to render it unobservable, as far as we consider mass loss rates of $10^{-8}~M_\odot$yr^-1, and FEBs identical to those assumed in the $\beta\:$Pic case, just scaling the evaporation parameters according to Fig. 4. These hypotheses need now to be discussed.

3.3 Lowering the mass loss rate

The value assumed in the simulations of Fig. 7 ( $10^{-8}\,M_{\odot}$yr^-1) falls in the range of values commonly assumed for Herbig Ae/Be stars. Of course, assuming higher values (such as $10^{-7}\,M_{\odot}$yr^-1, a still common value for such stars) leads to an even more drastic situation, the Mg II cloud becoming even thinner. Conversely, it is of valuable interest to find the maximum mass loss rate that would allow the FEBs to generate detectable variable components. This is illustrated in Figs. 10 and 11. Figure 10 is the same as Fig. 7b, but with $\dot{M}=10^{-11}M_\odot$yr^-1. Figure 11 displays only the Mg II k spectral signature of the FEB for various other $\dot{M}$values less than $\dot{M}=10^{-8}M_\odot$yr^-1. The lower the mass loss rate, the further away the wind shock surfaces are located from the FEB nucleus, and thus the deeper the variable components appear. Note that in all these runs, the production rate of the FEB was unchanged. The depth of the variable component appears in fact fully controlled by the wind.

Basically, for yr^-1, a detectable component is observable. The wind shock surfaces around the nucleus are located sufficiently far away from the FEB nucleus to allows the Mg II cloud to be large enough to have a significant spectral signature. The variable component turns out to be very deep for yr^-1, up to $\sim$75% of the continuum (Fig. 11). The size of the Mg II clouds is controlled by the wind, i.e. by the location of the tangential discontinuity, down to $\dot{M}\simeq5\,10^{-12}M_\odot$yr^-1; for weaker winds, the radiation pressure blows them away before they feel the wind, so that the spectral signature of the cloud does no longer grow. In this regime, we are back to a $\beta\:$Pic-like situation.

This $10^{-10}~M_\odot$yr^-1 wind threshold (valid for the FEB production rate assumed) is well below the usual range for Herbig Ae/Be stellar winds, so that we may conclude that the FEBs under consideration should not generate observable transient absorption events like those observed, at least in the environments of spherically symmetrical winds.

3.4 Increasing the production rate of the FEB

The depth of the variable feature generated by a passing FEB depends actually on the location of the wind shock surfaces around the nucleus. This location is controlled by the balance between the wind and the cometary mass fluxes, i.e. the production rate of the FEB (see Appendix). Hence lowering the wind at constant production rates leads to larger clouds, and consequently to deeper features. Conversely, increasing the production rates (dust and volatile) of the FEB should lead to a similar result, at least qualitatively. This is illustrated in Figs. 12 and 13. Figure 12 is equivalent to Fig. 7b, but the production rates have been multiplied by 1000 without changing the wind. We see that the wind shocks appear now further away from the nucleus, so that the Mg II cloud is thicker, leading to a very deep transient absorption. Note also that there is virtually no tail. Thanks to the high production rates, the Mg II cloud is (from a spectral point of view) optically thick; self-opacity causes thus most ions to feel virtually no radiation pressure, so that they remain around the nucleus for a longer time. This is also the reason why the resulting absorption component appears deeper than in Fig. 10 (the $\dot{M}=10^{-11}M_\odot$yr^-1 case), although the Mg II clouds achieve similar sizes in both cases.

	Figure 12: Same as Fig. 7b, but now with dust and volatiles production rates multiplied by 1000. A deep transient absorption component is generated
Open with DEXTER

	Figure 13: The spectral signature in the Mg II k line of the same FEB as in Figs. 7 and 10, for production rates multiplied by 10 (left) and 100 (right) with respect to the case of Fig. 7b
Open with DEXTER

Figure 13 shows the resulting Mg II k feature for intermediate cases, i.e. for production rates multiplied by 10 and 100 with respect to Fig. 7b. We note that a small feature is present in the $\times 10$case, but that it becomes really observable in the $\times 100$ case. This result holds for a wind mass loss rate $\dot{M}=10^{-8}~M_\odot$yr^-1. For $\dot{M}=10^{-7}~M_\odot$yr^-1, we must now multiply the production rates by 10000 to get observable features, and only by 10 for $\dot{M}=10^{-9}~M_\odot$yr^-1.

4 Discussion

The present simulations show that in the environment of a Herbig Ae/Be star with axisymmetrical wind, potential FEBs may generate observable spectral absorptions only i) if the stellar wind is less than $\sim$ $10^{-10}~M_\odot$yr^-1, or ii) if we increase the evaporation rates of the FEB by a factor 100 with respect to those scaled from the $\beta\:$Pic case (which were already 3 times higher than those deduced for $\beta\:$Pic). In order to get variable components as deep as in Fig. 1, the required factor is 500. An intermediate case is to consider $M= 10^{-9}~M_\odot$yr^-1. In that case, the production rates of the FEB have to be increased only by a factor 10 to get observable components, and 50 to get components as deep as in Fig. 1. We shall discuss now whether such situations are realistic or not.

4.1 Increasing the evaporation rate

In presence of a strong wind, the only way to get deep FEB spectral signatures is to significantly increase their evaporation rate. With such a high evaporation rate, the main problem concerns the survival of bodies down to the periastron values needed. In the case of $\beta\:$Pic, the dynamical characteristics of the FEBs are well modeled assuming they originate for the 4:1 mean-motion resonance with a Jovian-like planet (Beust & Morbidelli 1996, 2000).

This model appears to fit well the data for $\beta\:$Pic, but we want to point out that the proposed mechanism is extremely generic and should be active as soon as a planet on a moderately eccentric orbit ( ) is present, which is not in itself a strong requirement. The key point here for our purpose is that the process that brings the bodies from low-eccentricity to star-grazing orbits is gradual. In this scenario, a given body is able to become a star-grazer within $\sim$10⁴ revolutions of the perturbing planet, i.e. typically a few $10^5\,$yr depending on the orbital period of the planet. Hence if the dust begins to evaporate at $\sim$2AU from the star, then a body seen with a periastron of 0.04AU should have already undergone typically a few thousands of periastron passages within the dust evaporating zone before reaching this periastron value. It is not obvious whether it should survive to all these passages.

Figure 14: The minimum periastron that a body can reach before being fully evaporated, as a function of its radius, in the $\beta\:$Pic environment (thin line) and that of a B9 type star (thick lines), assuming evaporation rates 100 times larger in the latter case. The assumed dynamical evolution of the body is given by numerical integration of a particle trapped in the 4:1 resonance with a Jovian-like planet. In both cases, the rocky body is assumed to start evaporating when it reaches the dust evaporation limit

Open with DEXTER

In Beust & Morbidelli (2000), we investigated this question in the $\beta\:$Pic context, i.e. with a dust evaporating zone limited to $\sim$0.5AU and appropriate evaporation rates. The conclusion was that the FEBs need to be initially $\sim$15 km large (or more) in order to resist down to genuine star-grazing periastron values, as some observations (HVFs) seem to require. This of course does not exclude the presence of (numerous) smaller bodies, but they do not penetrate very deep into the dust evaporating zone, as they fully evaporate before. Now, in the Herbig Ae/Be case, the situation is even more drastic, as the dust evaporation zone extends 4 times further away from the star, and that the mere observation of variable features despite the presence of a strong wind seems to require the evaporation rates to be at least 100 times larger that the expected $\beta\:$Pic ones at the same stellar distances.

We apply here the same analysis to the B9 star case, i.e. with a dust evaporating limit of 1.8AU and evaporation rates 100 times larger as for $\beta\:$Pic. The compared result with $\beta\:$Pic is illustrated in Fig. 14, which is similar to Fig. 7 of Beust & Morbidelli (2000). Assuming a dynamical evolution given by numerical integration of a particle trapped in 4:1 mean-motion resonance with a perturbing Jovian-like planet (located at 10AU), we compute the mass loss of bodies of various initial radii at each periastron passage, as soon as they get into the dust evaporation zone. For various body radii, we thus compute the minimum periastron reached, which is is the periastron the body assumes when it turns out to be fully evaporated. Of course, larger bodies are expected to resist down to smaller periastron values, and the same body resists much less in the environment of a B9 star than around $\beta\:$Pic. This shows up in Fig. 14. As a matter of fact, we note from Fig. 14 that in the $\beta\:$Pic environment, bodies larger than $\sim$15 km are indeed able to get periastron values significantly less than the dust evaporation limit, while smaller bodies do not. In the B9 type star environment, this threshold turns out to be $\sim$100 km (Fig. 14). Therefore, for a FEB to be able to generate an absorption event when crossing the line of sight at 0.04AU from the star, it needs to be at least 100 km large, otherwise it is destroyed before. The constraint is even more drastic if we consider wind mass loss values in the range 10^-8- $10^{-7}~M_\odot$yr^-1.

In the peculiar case of HD 100546, if we assume similar winds, the evaporation rate of the FEB needs to be $\sim$500 times as large as for $\beta\:$Pic in order to reproduce the observed variable features (Fig. 1). This leads to a minimum size of $\sim$200 km instead of 100 km.

Assuming we observe towards Herbig Ae/Be FEBs that are at least $\sim$100 km-sized instead of $\sim$15 km in the case of $\beta\:$Pic may help understanding the discrepancy concerning the evaporation rates. The ${\rm d}E/{\rm d}t$ values plotted in Fig. 4 are the evaporation rates per unit surface, i.e. they compare the global evaporation rates between the various stellar environments for identical bodies. Between a B9 star and $\beta\:$Pic, the ratio is 3. Now, if the typical Herbig Ae/Be FEB is 100 km-sized, while the one for $\beta\:$Pic is 15 km, the global evaporation rate ratio between the two bodies should by $(100/15)^2\times3\simeq 44$, which is comparable to what we have deduced from our simulations.

The problem with this scenario is the number of such bodies, and consequently the amount of mass in the disk that should be necessary in order to maintain a statistically detectable spectral FEB activity. Note also that if there are many 100 km-sized bodies, there should also be many more smaller ones. We investigated the relevance of this question to the $\beta\:$Pic case in Beust & Morbidelli (2000), coming to the conclusion that the total mass of planetesimals necessary in the disk for sustaining the present FEB activity (several hundreds of events per year) while refilling the 4:1 resonance by mutual collisions among planetesimals or planet migration was at least $10\,M_\oplus\,$ per AU . This is only marginally compatible with realistic estimates. Obviously a more refined study of this question is necessary.

In the present case of Herbig Ae/Be stars, it is almost impossible to derive any comparable reliable estimate, as we still do not know the spectral events frequency. All we can say is that the events are not rare, as they have been observed. In the case of HD 100546, we can point out that variable events have been observed more than once (Grady et al. 1997; Viera et al. 1999). On the basis of observed variation time-scales of a few tens of hours, Grady et al. (2000) give nevertheless tentative events frequency estimates for some stars else than HD 100546, with results ranging between one event per day and one every five days. This is $\sim$10 times less than the estimates for $\beta\:$Pic. In Beust & Morbidelli (2000), we inferred from dynamical simulations a population of bodies $N\simeq10^8$-10⁹ bodies per AU in the vicinity of the 4:1 resonance. The mass density $M_{{\rm l}}$ within the disk then reads

$$M_{{\rm l}} = \frac{20}{3}\pi Nr_{{\rm min}}^3\rho\left( \sqrt{\frac{r_{{\rm max}}}{r_{{\rm min}}}}-1\right),$$ (13)

where $r_{{\rm min}}$ and $r_{{\rm max}}$ are the minimum and maximum radii for candidate FEBs, and rho their density. Assuming for $\beta\:$Pic $r_{{\rm max}}=500\,$km, $r_{{\rm min}}=15\,$km, and $\rho=1\,$gcm^-3, we derive values for $M_{{\rm l}}$ ranging between 6.8 and 68$M_\oplus$ per AU. Transposing now this estimate to the Herbig Ae/Be case, i.e. with $r_{{\rm min}}=100\,$km and N 10 times less than for $\beta\:$Pic, we derive

$$520\,M_\oplus\,\mbox{AU}^{-1}\mathrel{\mathchoice {\vcenter{\... ...yle ...$$ (14)

which is now unrealistic. The estimate turns out to be even higher if we assume $r_{{\rm min}}=200\,$km as it seems to be required for explaining deep events like displayed in Fig. 1.

In the case of AB Aur, Grady et al. (1999) stress that the infall activity is even higher than for $\beta\:$Pic and HD 100546, arguing that signatures of infalling gas was detected in all high-resolution UV spectra of AB Aur from 1978 to 1996. Consequently, the disk population and mass estimates will be even higher than for HD 100546, which turns out to be even more unrealistic.

We may then stress that even if we cannot derive any precise events frequency estimate, if we assume that these events are due to 100 km-sized bodies, the inferred "planetesimal'' disk mass turns out to be truly unrealistic.

4.2 Lowering the mass loss rate: Specific stellar cases

If we want to get detectable FEB features while keeping reasonable evaporation rates, we need to lower the mass loss rate of the star. This issue needs now to be discussed for each stellar case.

In the case of AB Aur, where the wind has been successfully modeled as spherically symmetric with a mass loss rate of $1.8\,10^{-8}~M_\odot$yr^-1, the conclusion is straightforward: FEBs, if present, could probably not be observed in absorption lines. This conclusion holds in fact for all stars where similar mass loss rates have been measured, such as those modeled together with AB Aur by Bouret & Catala (1998). It can in fact be extended to all Herbig Ae/Be stars of comparable ages (a few 10⁶yr), as they are all expected to have strong winds. This conclusion is independently supported by recent observation of an accretion episode towards UX Ori by Natta et al. (2000), showing that the infalling gas is compatible with solar rather than cometary (i.e., volatile depleted) composition. For such stars, this suggests to attribute the accretion episodes observed to magnetically driven accretion columns (see e.g. models by Hatman et al. 1994) rather than to FEBs. For AB Aur however, we mentioned above that the wind might not actually be spherically symmetric, but only modeled as spherically symmetric thanks to its high viewing latitude from the Earth. But even if this is true, in that case this means that we are actually viewing the star across the wind, so that the main conclusion concerning the detectability of FEB signatures still holds.

This conclusion does not straightforwardly extend to older Herbig Ae/Be stars ($\sim$10⁷yr) such as HD 100546. For HD 100546 specifically, Viera et al. (1999) conclude on the basis of H$_\alpha$ profiles to the presence of at least three distinct components in the circumstellar envelope of the star: a discrete, variable accretion close to the star, a more stable one in the remote envelope, and a wind that appeared stable during the period they observed the star. They however do not derive any mass loss rate. In fact, modeling this structure appears here more complex than for P-Cygni class Herbig Ae/Be stars like AB Aur, for which the wind component dominates the line profile.

On this simple basis, we may stress that the mass loss rate of HD 100546 should not exceed that of AB Aur. Based on a statistical study of T Tauri stars, Muzerolle et al. (2000) showed recently that accretion rates drop by 2 orders of magnitude between ages of 10⁶ and $10^7\,$yr. This holds for T Tauri stars, and for accretion. If we were to believe this to be still true for Herbig Ae/Be stars, and to similarly apply for wind mass loss rates, then scaling from measured values for young Herbig Ae/Be stars like AB Aur, we would derive a mass loss rate estimate for HD 100546 of a few $10^{-10}~M_\odot$yr^-1. This is close to the limit under which FEB signatures begin to be observable. We may consider this value as a lower limit. Therefore, provided this lower limit applies, FEBs could be marginally observable around stars like HD 100546, although strong components like those of Fig. 1 are hardly reproduced with $10^{-10}~M_\odot$yr^-1.

Up to now, we have considered a spherically symmetric wind geometry. If this approximation may be considered as valid for young Herbig Ae/Be stars like AB Aur, it is far from being obvious to apply for older stars like HD 100546. Viera et al. (1999) suggest conversely that the regions of the stellar environment of HD 100546 with wind activity should be limited to the higher latitudes, while lower latitude regions should be wind free, accretion of matter taking place close to the equatorial plane. In some specific cases such as HD163296, there is observational evidence for a collimated outflow, i.e., far from spherically symmetric (Devine et al. 2000).

Such non-spherically symmetric geometries can be generated magnetically. This is especially the case if the stars holds a dipole-like fossil magnetic field. In that case, models by Paatz & Camenzind (1996) and Strafella et al. (1998) show that in the first few stellar radii above the stellar surface, the dipole-like magnetic field is expected to co-rotate with the star, leaving a wind free cavity at low latitude, while at higher latitudes the lines should be open and thus transport matter away.

If this picture is correct for HD 100546, then we should expect to potentially observe spectral signatures of FEBs if they cross the line of sight inside the wind free cavity, i.e., at low latitude. As a matter of fact, the mean-motion resonance model invoked at the source for the FEB phenomenon towards $\beta\:$Pic predicts that most of the FEBs should remain at moderate orbital inclination with respect to the equatorial plane of the star.

Besides, we know now that HD 100546 is not viewed edge-on. Pantin & Lagage (2001) and Augereau et al. (2001) report recently and independently the detection of a dusk disk surrounding the star, in thermal emission for Pantin & Lagage (2001) and in scattered light for Augereau et al.(2001). Both detections agree for a disk inclination of $40\hbox{$^\circ$ }$ with respect to pole-on, based on the shape of the circumstellar nebulosity. This means that the line of sight is located at $\sim$ $50\hbox{$^\circ$ }$ stellar latitude. It is thus far from being obvious to state whether at $50\hbox{$^\circ$ }$ latitude we are viewing the stellar surface across the wind or not. Moreover, even if the wind is not strong enough to prevent FEB detection, it is still not obvious whether FEBs at $\sim$ $50\hbox{$^\circ$ }$ inclination can be dynamically generated. This specific issue will be addressed in forthcoming work.

In this framework, even a mass loss rate measurement does not allow to state whether FEBs can be observed or not towards such a star, as this depends on the geometry of the wind. Specific, more sophisticated modeling of this stellar environment is therefore required. Moreover, a more regular survey of the spectrum of stars like HD 100546 would help to refine the statistics on the occurrence of these FEB-like events.

5 Conclusions

Contrary to the $\beta\:$Pic case, interpreting the transient spectral events observed in metallic lines of several Herbig Ae/Be stars in the frame of the Falling Evaporating Bodies scenario appears much less straightforward, as the presence of strong stellar winds around such stars makes it much more difficult to actually observe deep absorption components generated by individual FEBs. Besides, in the environment on such stars, the refractory material start evaporating much further away from the star than around $\beta\:$Pic, which can increase the probability of FEB detection. However, the observations showing high velocity absorption features are compatible with FEBs crossing the line of sight at a less than 0.1AU only. In such conditions, and assuming that the process that brings them down to such small stellar distances is the same as the one invoked for $\beta\:$Pic FEBs, the only way make them i) resist a sufficiently long time and ii) to render their spectral signature easily detectable in presence of a strong spherically symmetric wind typical for Herbig Ae/Be stars is to assume that they are at least $\sim$100 km-sized instead of 15 km for $\beta\:$Pic FEBs, and that their evaporation rate should be increased with respect to $\beta\:$Pic FEBs proportionally to the flux of the star and to their surface.

This assumption immediately raises the question whether the amount of mass that should be involved in this process is realistic or not. It is still impossible to definitely answer this question, because of the lack of observational statistics on each of these stars. The first events frequency estimates on some of them lead nevertheless to unrealistic masses. A more detailed analysis and a better observational coverage will help solving this question.

It is nevertheless possible to stress that for young Herbig Ae/Be stars like AB Aur with well identified spherically symmetric winds with mass loss rates of $\sim$10 $^{-8}~M_\odot\,$yr^-1 or above, the FEB hypothesis is probably not realistic to account for variable accretion events.

For older Herbig Ae/Be stars like HD 100546 however, the conclusion is not so straightforward, as winds could be close to the limit under which FEB signatures begin to be observable, and because various observations and modeling suggest that the wind is confined to higher stellar latitudes. Hence FEBs could still be observed despite a strong stellar wind if they cross a the line of sight at low stellar latitude. Investigating this issue will require i) more precise measurements of mass loss rates, ii) more extensive modeling adapted to each stellar case, and iii) a more regular survey of the variable spectral events occurring towards such stars.

Appendix A Stellar wind-comet interaction formalism

The formalism we assume for the interaction between the stellar wind and the cometary plasma is classical and simple, mainly taken from Galeev et al. (1985) and Gombosi et al. (1983). In the referential frame of the comet, the wind is assumed to come from infinity, with a plan-parallel geometry. We first describe the flow along the axis cometary nucleus - star. The continuity, momentum, and energy conservation equations for the mass loaded solar wind along this line, in front of the bow shock, read

   

$$\frac{{\rm d}}{{\rm d}r}(\rho u) -\frac{mn}{\tau}\;;$$ = (A.1)

$$\frac{{\rm d}}{{\rm d}r}\left(\rho u^2+P+\frac{B^2}{2\mu_0}\right) 0\;;$$ = (A.2)

$$\frac{{\rm d}}{{\rm d}r}\left[u\left(\frac{\rho u^2}{2}+\frac{P\gamma}{\gamma-1} +\frac{B^2}{2\mu_0}\right)\right] 0\;.$$ = (A.3)

In these equations, u stands for the velocity of the wind with respect to the comet, n for the number density of the cometary neutrals, $\tau$ for their ionization time-scale, m for their mass; $\rho$ is the mass density of the wind, r is the distance to the cometary nucleus, B is the magnetic field strength and $\gamma$ is the adiabatic index of the gas. The number density of cometary neutrals, in a spherical expansion, reads

$$n=\frac{Q}{4\pi r^2 v_{{\rm e}}}(1-\xi(r))\;,$$ (A.4)

where Q is the production rate of the comet, $v_{{\rm e}}$the expansion velocity of the cometary neutrals, and $\xi$ is the ionization level of the volatile gas. If $\xi$ remains small, it may be taken as zero in Eq. (A.4), and subsequently in Eq. (A.1). This case is often referred as the linear comet case, which holds close enough to the nucleus but not at larger distances. If we take into account the decrease of the neutral density due to photoionization, a more correct and frequently used expression of the neutral density is

$$n=\frac{Q}{4\pi r^2 v_{{\rm e}}}\exp\left(-\frac{r} {v_{{\rm e}}\tau}\right).$$ (A.5)

However, the latter expression is valid only if we may neglect the recombination of ions. As pointed out by Galeev et al. (1985), this is certainly valid in the coma of Solar comets, as the density is fairly low. The FEBs we are considering have outgassing rates at least 10⁴ times higher than typical Solar comets, so that we think that we may not neglect recombination. A correct description is given in Gombosi et al. (1983). In a spherically symmetric expansion model at constant velocity $v_{{\rm e}}$, the differential equations governing the densities of neutrals n(r) and ions $n_{{\rm i}}(r)$are

$$\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2v_{{\rm e}}n(r)\right) -\frac{n(r)}{\tau}+\alpha n_{{\rm i}}^2(r)\;,$$ = (A.6)

$$\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2v_{{\rm e}}n_{{\rm i}}(r)\right) \frac{n_{{\rm i}}(r)}{\tau}-\alpha n_{{\rm i}}^2(r)\;,$$ = (A.7)

where $\alpha$ is the recombination rate coefficient, assumed constant here. Combining with $n_i(r)+n(r)=Q/(4\pi r^2v_{{\rm e}})$, we derive the differential equation for the ionization level $\xi(r)=n_{{\rm i}}(r)/(n(r)+n_{{\rm i}}(r))$:

$$\frac{{\rm d}\xi(r)}{{\rm d}r}=\frac{1}{v_{{\rm e}}\tau}-\fra... ...m e}}\tau}-\frac{\alpha Q\xi^2(r)}{4\pi r^2v_{{\rm e}}^2}\cdot$$ (A.8)

If we introduce now the dimensionless variable and parameter

$$s=\frac{r}{v_{{\rm e}}\tau}\qquad\mbox{and} \qquad G=\frac{\alpha Q}{4\pi v_{{\rm e}}^3\tau},$$ (A.9)

this equation reduces to

$$\frac{{\rm d}\xi(s)}{{\rm d}s}=1-\xi(s)-G\frac{\xi^2(s)}{s^2}\cdot$$ (A.10)

G should be taken as non-zero only where ions are actually present for recombining, i.e., inside the tangential discontinuity. Outside, each newly produced ion is immediately picked-up by the stellar wind, so that it cannot take part to recombination. Hence we must assume G=0 as soon as we are considering the flow inside the wind (i.e., outside the ionopause), and $G\neq 0$ where it is protected from the wind.

If we set G=0 (no recombination), the solution of this equation is $\xi(s)=1-{\rm Cst.}\times{\rm e}^{-s}$. If we did not consider the zone inside the tangential discontinuity where recombination is able to take place, we would add the boundary condition $\xi(0)=0$ (we neglect the physical radius of the nucleus), i.e. ${\rm Cst.}=1$, add this would bring us back to Eq. (A.5). Here the constant will be fixed by continuity with the solution inside the tangential discontinuity at the stagnation point of the wind $r=R_{{\rm stag}}$ (or $s=S_{{\rm stag}}$).

Inside the ionopause, where $G\neq 0$, the solution of Eq. (A.10) is more complex. As shown by Gombosi et al. (1983), the general solution to this equation involves the Whittaker function, but with the boundary condition $\xi(0)=0$, the solution may be expressed as

$$\xi(s)=\frac{2s\,I_{K}(s/2)}{(2K+s+1)\,I_{K}(s/2)+s\,I_{K+1}(s/2)},$$ (A.11)

where $I_\nu$ is the modified Bessel function of the first kind and order $\nu$, and where $K=\sqrt{1+4G}/2$. Now, a fairly accurate fit to this function is

$$\xi(s)=\frac{2s}{\sqrt{P^2+4s^2}}-\frac{2Ps^2(P-4)}{(P^2+4s^2)^2},$$ (A.12)

with P=2K+1. This fit has the correct asymptotic behaviors at s=0 and s=1, and apart for G values less that unity (which shall never be the case in our applications), the maximum relative error between this fit and the actual solution never exceeds 6%. Hence we decide to use the fit (A.12) in all our applications. Following Gombosi et al. (1983), we also decided to assume $\alpha=10^{-6}\,{\rm cm}^3\,{\rm s}^{-1}$ in the numerical applications.

The solution outside the ionopause then reads by continuity

$$\xi(s)=1-\left(1-\xi\left(S_{{\rm stag}}\right)\right) \exp\left(S_{{\rm stag}}-s\right),$$ (A.13)

where $\xi\left(S_{{\rm stag}}\right)$ is obtained from Eq. (A.12).

The continuity Eq. (A.1) can be now integrated. For linear comets first, we find

$$\hat{\rho}\hat{u}\equiv \left(\frac{\rho}{\rho_\infty}\right)... ...Qm}{4\pi v_{{\rm e}}\tau\rho_{\infty}u_{\infty}}\frac{1}{r}\;,$$ (A.14)

where $\infty$-subscribed quantities denote the corresponding quantities at infinity, i.e. in the unperturbed wind. $\hat{\rho}$and $\hat{u}$ are the same parameters normalized to their corresponding unperturbed values (see Galeev et al. 1985; Flammer 1991). If we take into account the ionization level taken from Eq. (A.13), the result is

$$\hat{\rho}\hat{u}= 1-\frac{Qm\left(1-\xi\left(S_{{\rm stag}}\... ...-\!\int_{s}^{+\infty} \frac{{\rm e}^{-u}}{u}\,{\rm d}u\right].$$ (A.15)

If we assume that in the supersonic domain, the flow remains super-Alfvénic, then the magnetic terms in Eqs. (A.2) and (A.3) can be neglected, and we can solve these equations to derive $\hat{\rho}$, $\hat{p}$ and $\hat{u}$ as a function of $\hat{\rho}\hat{u}$, itself given by Eq. (A.15):

 

$$\hat{u} \frac{1}{\hat{\rho}\hat{u}}\frac{\gamma}{\gamma+1} \left(1\pm\sqrt{1-\frac{\gamma^2-1}{\gamma^2}\hat{\rho}\hat{u}}\right)\;,$$ =

$$\hat{\rho} \frac{\hat{\rho}\hat{u}}{\hat{u}},\;\hat{p}\;=\;1- (\hat{\rho}\hat{u})\hat{u}\;.$$ = (A.16)

In this expression, the plus sign in front if the square root should apply to a supersonic flow, and the minus sign to a subsonic flow. This results shows that this flow regime can only be possible as long as $\hat{\rho}\hat{u}\leq\gamma^2/(\gamma^2-1)$. If the mass loading is important enough (this is usually the case), then a shock must form to decelerate the wind down to subsonic. Dedicated simulations Schmidt & Wegmann (1982) showed that this shock was expected to occur at Mach number $M\simeq 2$. This is the value we will assume. Following Galeev et al. (1985), we will also assume $\gamma=2$ in the supersonic flow. This in turn allows to calculate the normalized mass flux at the shock $(\hat{\rho}\hat{u})_{{\rm s}}$:

$$(\hat{\rho}\hat{u})_{{\rm s}}=\frac{\gamma^2}{\gamma-1} \frac{M^2\left[M^2(\gamma^2-1)+1\right]}{(M^2\gamma+1)^2}\;.$$ (A.17)

The plasma parameters $(u_2,\rho_2,p_2)$ just behind the shock are then computed assuming $\gamma=5/3$ and selecting the minus sign in Eq. (A.16), with $(\hat{\rho}\hat{u})_{{\rm s}}$ (see Galeev et al. 1985 for details). After the shock, the flow is expected to undergo a rapid isotropization of the velocity distribution, so that Eqs. (A.1)-(A.3) no longer apply. Equations (A.1) and (A.2) must be replaced, in the linear comet case (this is a valid approximation well inside the coma), by

  

$$\frac{1}{A}\frac{{\rm d}}{{\rm d}r}(\rho uA) -\frac{Qm}{4\pi v_{{\rm e}}\tau r^2}\;;$$ = (A.18)

$$\frac{1}{A}\frac{{\rm d}}{{\rm d}r}(\rho u^2A)+\frac{{\rm d}p}{{\rm d}r} -\frac{{\rm d}}{{\rm d}r}\left(\frac{B^2}{2\mu_0}\right)- \frac{B^2}{\mu_0r_{{\rm c}}}\;,$$ = (A.19)

where A is the flow-tube area and $r_{{\rm c}}$ is the radius of curvature of the magnetic filed lines. In the limiting case of negligible magnetic effects and gas cooling, and assuming $uA={\rm cst.}$along the flow, Galeev et al. (1985) show that Eq. (A.19) leads to the very simple solution $\rho(r)u(r)^{3/4}={\rm cst.}$ along the flow. Combined with Eq. (A.19) we derive the velocity profile

$$u(r)=u_2\left[1-\frac{R_{{\rm L}}}{r} \left(1-\frac{r}{R_{{\rm s}}}\right)\right]^4\;,$$ (A.20)

where

$$R_{{\rm L}}=\frac{Qm}{12\pi v_{{\rm e}}\tau\rho_2 u_2}\cdot$$ (A.21)

The flow velocity then decreases as r decreases, down to the stagnation point $r=R_{{\rm stag}}$ where u(r) vanishes. We thus get

$$\frac{1}{R_{{\rm stag}}}=\frac{1}{R_{{\rm L}}}+\frac{1}{R_{{\rm s}}}\,$$ (A.22)

which enables to derive the location of the tangential discontinuity $R_{{\rm stag}}$ as a function the location of the bow shock $R_{{\rm s}}$. Now setting $\hat{\rho}\hat{u}=(\hat{\rho}\hat{u})_{{\rm s}}$ in Eq. (A.15), we are able to numerically solve this equation to derive $R_{{\rm s}}$, where $R_{{\rm stag}}$ is obtained from Eq. (A.22) and the ionization level at this point from Eq. (A.12).

This description holds for the flow along the comet - star axis. Now we want to describe the flow off the axis. We first assume an axial symmetry of the flow around that axis, so that any planar description is valid. We choose to describe the flow in a plane containing the comet - star axis, and we define a referential frame (OXY) centered on the nucleus where the OX axis points towards the star and where OY is perpendicular to that direction. We will also neglect for simplicity any motion of the gas that is not parallel to OX, assuming in fact that the radial motion dominates the transversal one. We will then describe the flow along an axis parallel to OX, passing at distance yfrom the nucleus. In the supersonic part of the flow, Eqs. (A.1)-(A.3) still apply, provided ${\rm d}/{\rm d}r$ must be replaced by ${\rm d}/{\rm d}x$. The density of neutrals now reads for non-linear comets:

$$n=\frac{Q}{4\pi (x^2+y^2)v_{{\rm e}}} \left(1-\xi\left(\sqrt{x^2+y^2}\right)\right).$$ (A.23)

If we restrict ourself to the limiting case of linear comets, i.e., setting $\xi\equiv0$(otherwize the equation cannot be integrated easily), the revised Eq. (A.2) can be solved to yield

$$\hat{\rho}\hat{u} =1+\frac{Qm}{4\pi v_{{\rm e}}\tau\rho_{\infty}u_{\infty}} \frac{1}{y}\arctan\left(\frac{y}{x}\right).$$ (A.24)

The solution of Eqs. (A.2) and (A.3) remains unchanged, like the normalized mass flux on the shock $(\hat{\rho}\hat{u})_{{\rm s}}$. In the linear comet formalism, Eq. (A.24) may be rewritten as

$$\hat{\rho}\hat{u}=1+R_{{\rm s}}\left[(\hat{\rho}\hat{u})_{{\rm s}}-1 \right]\frac{1}{y}\arctan\left(\frac{y}{x}\right),$$ (A.25)

where $R_{{\rm s}}$ was introduced above. Although this strictly holds for linear comets, we decided to use it for our non-linear comets to describe the flow, provided $R_{{\rm s}}$ is derived solving Eq. (A.15). The location of the bow shock for every y is defined by $\hat{\rho}\hat{u}=(\hat{\rho}\hat{u})_{{\rm s}}$. Introducing this into Eq. (A.25) gives the Cartesian planar equation of the shock:

$$\frac{1}{R_{{\rm s}}}=\frac{1}{y}\arctan\left(\frac{y}{x}\right).$$ (A.26)

This is better rewritten in polar coordinates $(x=r\cos\theta,y= r\sin\theta)$ as

$$r=R_{{\rm s}}\frac{\theta}{\sin\theta}\;,$$ (A.27)

showing that the shock assumes as expected a parabolic-like shape. If we turn now to describe the flow after the shock, we need to rewrite Eqs. (A.18) and (A.19). As for the supersonic flow, the momentum Eq. (A.19) remains unchanged, once ${\rm d}/{\rm d}r$ has been changed by ${\rm d}/{\rm d}x$; the same applies to the continuity Eq. (A.18), apart from the fact that we need to replace r² by x²+y². As the momentum equation does not change, we still derive $\rho(x)u(x)^{3/4}={\rm cst}$ along the flow. Combined with the revised continuity equation, Eq. (A.20) becomes now

$$u(x)=u_2\left[1-\frac{R_{{\rm L}}}{x}\left( \frac{x}{y}\arctan\left(\frac{y}{x}\right)-\frac{x}{R_{{\rm s}}}\right) \right]^4.$$ (A.28)

The tangential discontinuity is again defined as the locus of the points where u(x) vanishes, which turns out to verify

$$\frac{1}{R_{{\rm stag}}}=\frac{1}{y}\arctan\left(\frac{y}{x}\... ...\mbox{or}\qquad r=R_{{\rm stag}}\frac{\theta}{\sin\theta}\cdot$$ (A.29)

Hence the tangential discontinuity assumes the same shape as the bow shock, but is it much narrower, as $R_{{\rm stag}}\ll R_{{\rm s}}$. The shapes of these surfaces is represented in Fig. 6.

References

 
• Artymowicz, P. 1988, ApJ, 335, L79 In the text NASA ADS  
• Artymowicz, P. 1997, Ann. Rev. Earth Planet. Sci., 25, 175 In the text  
• Augereau, J.-C., Lagrange, A.-M., Mouillet, D., & Ménard F. 2001, A&A, submitted In the text  
• Babel, J. 1995, A&A, 301, 823 In the text NASA ADS  
• Bertout, C. 1989, ARA&A, 27, 351 In the text NASA ADS  
• Beust, H., Lagrange-Henri, A.-M., Vidal-Madjar, A., & Ferlet, R. 1989, A&A, 223, 304 In the text NASA ADS  
• Beust, H., Lagrange-Henri, A.-M., Vidal-Madjar, A., & Ferlet, R. 1990, A&A, 236, 202 In the text NASA ADS  
• Beust, H., & Morbidelli, A. 1996, Icarus, 120, 358 In the text NASA ADS  
• Beust, H., Lagrange, A.-M., Plazy, F., & Mouillet, D. 1996, A&A, 310, 181 In the text  
• Beust, H., Lagrange, A.-M., Crawford, I. A., et al. 1998, A&A, 338, 1015 In the text NASA ADS  
• Beust, H., & Morbidelli, A. 2000, Icarus, 143, 170 In the text NASA ADS  
• Boggess, A., Bruhweiler, F. C., Grady, C. A., et al. 1991, ApJ, 377, L49 In the text NASA ADS  
• Bouret, J.-C., Catala, C., & Simon, T. 1997, A&A, 328, 606 In the text NASA ADS  
• Bouret, J.-C., & Catala, C. 1998, A&A, 340, 163 In the text NASA ADS  
• Bouret, J.-C., & Catala, C. 2000, A&A, 359, 1011 In the text NASA ADS  
• Bouwman, J., de Koter, A., van den Ancker, M. E., & Waters, L. B. F. M. 2000, A&A, 360, 213 In the text NASA ADS  
• Castor, J. I., Abott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 In the text NASA ADS  
• Catala, C., & Kunasz, P. B. 1987, A&A, 174, 158 In the text NASA ADS  
• Catala, C. 1989, Herbig, Ae and Be stars, in ed. B. Reipurth, Low-Mass Star Formation and Pre-Main-Sequence Objects, ESO Garching, 471 In the text  
• Clarke, D., Smith, R. A., & Yudin, R. V. 1999, A&A, 347, 590 In the text NASA ADS  
• Crifo, F., Vidal-Madjar, A., Lallement, R., Ferlet, R., & Gerbaldi, M. 1997, A&A, 320, L29 In the text NASA ADS  
• Damiani, F., Micela, G., Sciortino, S., & Harnden, F. R. 1994, ApJ, 436, 807 In the text NASA ADS  
• Devine, D., Grady, C. A., Kimble, R. A., et al. 2000, ApJ, 542, L115 In the text NASA ADS  
• de Winter, D., Grady, C. A., van den Ancker, M. E., Pérez, M. R., & Eiroa, C. 1999, A&A, 343, 137 In the text NASA ADS  
• Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 In the text NASA ADS  
• Flammer, K. R. 1991, The Global Interaction of Comets with the Solar Wind, ed. R. L. Newburn, J. Rahe, & M. Neugebauer, in Comets in the Post-Halley Era., vol. 2 (Kluwer Academic Publishers), 1125 In the text  
• Flammer, K. R., Mendis, D. A., & Shapiro, V. D. 1997, ApJ, 482, 1021 In the text NASA ADS  
• Galeev, A. A., Cravens, T. E., & Gombosi, T. I. 1985, ApJ, 289, 807 In the text NASA ADS  
• Gombosi, T. I., Horanyi, M., Kecskemety, K., Cravens, T. E., & Nagy, A. F. 1983, ApJ, 268, 889 In the text NASA ADS  
• Grady, C. A., & Silvis, J. M. S. 1993, ApJ, 402, L61 In the text NASA ADS  
• Grady, C. A., Pérez, M. R., & Thé, P. S. 1993, A&A, 274, L847 In the text  
• Grady, C. A, Pérez, M. R., Talavera, A., et al. 1996, A&AS, 120, 157 In the text NASA ADS  
• Grady, C. A., Pérez, M. R., Talavera, A., et al. 1996, ApJ, 471, L49 In the text NASA ADS  
• Grady, C. A., Sitko, M. L., Bjorkman, K. S., et al. 1997, ApJ, 483, 449 In the text NASA ADS  
• Grady, C. A., Pérez, M. R., Bjorkman, K. S., & Massa, D. 1999, ApJ, 511, 925 In the text NASA ADS  
• Grady, C. A., Woodgate, B., Bruhweiler, F. C., et al. 1999, ApJ, 523, L151 In the text NASA ADS  
• Grady, C. A., Sitko, M. L., Russell, R. W., et al. 2000, Infalling Planetesimals in Pre-Main Sequence Stellar Systems, in ed. V. Mannings, A. P. Boss, & S. S. Russell, Protostars and Planets IV (University of Arizona Press, Tucson), 613 In the text  
• Grinin, V. P., Thé, P. S., de Winter, D., et al. 1994, A&A, 292, 165 In the text NASA ADS  
• Grinin, V., Natta, A., Tambovtseva, L. 1996, A&A, 313, 857 In the text NASA ADS  
• Grinin, V. P., Kozlova, O. V., Thé, P. S., & Rostopchina, A. N. 1996, A&A, 309, 474 In the text NASA ADS  
• Hartmann, L., Hewett, R., & Calvet, N. 1994, ApJ, 426, 669 In the text NASA ADS  
• Hobbs, L. M. 1986, ApJ, 308, 854 In the text NASA ADS  
• Huebner, W. F., Boice, D. C., Schmidt, H. U., & Wegmann, R. 1991, Structure of the coma: Chemistry and Solar wind interaction, in Comets in the Post-Halley Era., vol. 2, ed. R. L. Newburn, J. Rahe, & M. Neugebauer (Kluwer Academic Publishers), 907 In the text  
• Krankowsky, D. 1991, The composition of comets, in Comets in the Post-Halley Era., vol. 2, ed. R. L. Newburn, J. Rahe, & M. Neugebauer (Kluwer Academic Publishers), 855 In the text  
• Kurucz, R. L. 1979, ApJS, 40, 1 In the text NASA ADS  
• Lagrange-Henri, A.-M., Ferlet, R., Vidal-Madjar, A., et al., 1990, A&AS, 85, 1089 In the text NASA ADS  
• Lagrange-Henri, A.-M., Beust, H., Ferlet, R., Hobbs, L. M., & Vidal-Madjar, A., A&A, 227, L13 In the text  
• Lagrange, A.-M., Plazy, F., Beust, H., et al. 1996, A&A, 310, 547 In the text NASA ADS  
• Lagrange, A.-M., Beust, H., Mouillet, D., et al. 1998, A&A, 330, 1091 In the text NASA ADS  
• Lagrange, A.-M., Backman, & D. E., Artymowicz, P. 2000, Planetary Material around Main Sequence Stars, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, S. S. & Russell (University of Arizona Press, Tucson), 639 In the text  
• Lamy, P. L. 1974, A&A, 35, 197 In the text NASA ADS  
• Laor, A., & Draine, B. T. 1993, ApJ, 402, 441 In the text NASA ADS  
• Lecavelier des Etangs, A., Vidal-Madjar, A., & Ferlet, R. 1996, A&A, 307, 542 In the text NASA ADS  
• Lecavelier des Etangs, A., Deleuil, M., Vidal-Madjar, A., et al. 1997, A&A, 325, 228 In the text NASA ADS  
• Malfait, K., Waelkens, C., Waters, L. B. F. M., et al. 1998, A&A, 332, L25 In the text NASA ADS  
• Mouillet, D., Larwood, J. D., Papaloizou, J. C. B., & Lagrange, A.-M. 1997, MNRAS, 292, 896 In the text NASA ADS  
• Muzerolle, J., Calvet, N., Briceño, C., Hartmann, L., & Hillenbrand, L. 2000, ApJ, 535, L47 In the text NASA ADS  
• Natta, A., Grinin, V. P., & Tambovtseva, L. V. 2000, ApJ, 542, 421 In the text NASA ADS  
• Paatz, G., & Camenzind, M. 1996, A&A, 308, 77 In the text NASA ADS  
• Pantin, E., & Lagage, P.-O. 2001, A&A, submitted In the text  
• Parker, E. 1958, ApJ, 128, 664 In the text NASA ADS  
• Rickman, H. 1991, Physico-Dynamical Evolution of aging Comets, in Interrelations between physics and dynamics for minor bodies in the Solar System, ed. D. Benest, & C. Froeschlé (Editions Frontières), 197 In the text  
• Schmidt, H. U., & Wegmann, R. 1982, in Comets, ed. L. L. Wilkening (Univ. Arizona Press, Tucson) In the text  
• Sitko, M. L., Grady, C. A., Lynch, D. K., Russell, R. W., & Hanner, M. S. 1999, ApJ, 510, 408 In the text NASA ADS  
• Sorelli, C., Grinin, V. P., & Natta, A. 1996, A&A, 309, 155 In the text NASA ADS  
• Smith, B., & Terrile, R. 1984, Sci, 226, 1421 In the text NASA ADS  
• Strafella, F., Pezzuto, S., & Corciulo, G. G. 1998, ApJ, 505, 299 In the text NASA ADS  
• van den Ancker, M. E., de Winter, D., & Tjin A Djie, H. R. E. 1998, A&A, 330, 145 In the text NASA ADS  
• van den Ancker, M. E., Bouwman, J., Wesselius, P. R., et al. 2000, A&A, 357, 325 In the text NASA ADS  
• Vidal-Madjar, A., Lagrange-Henri, A.-M., Feldman, P. D., et al. 1994, A&A, 290, 245 In the text NASA ADS  
• Vidal-Madjar, A., Lecavelier des Etangs, A., & Ferlet, R. 1998, Planet. Space Sci., 46, 629 In the text  
• Viera, S. L. A., Pogodin, M. A., & Franco, G. A. P. 1999, A&A, 345, 559 In the text NASA ADS  
• Waelkens, C., Waters, L. B. F. M., de Graauw, M. S., et al. 1996, A&A, 315, L245 In the text NASA ADS  
• Welsh, B. Y., Craig, N., Crawford, I. A., & Price, R. J. 1998, A&A, 338, 674 In the text NASA ADS