format: text

text:
# Long-term X-ray variability of the microquasar system LS 5039/RX J1826.2-1450

## 1 Introduction

LS~5039/RX~J1826.2$-$1450 was first identified as a new massive X-ray binary by Motch et al. (motch97 (1997)). Interest in this source has grown significantly because it has turned out to be a source of persistent radio-emitting relativistic jets, and because it is believed to be physically associated with the $`\gamma `$-ray source 3EG~J1824$-$1514 (Paredes et al. paredes00 (2000), paredes02 (2002)).

The first radio detection was reported by Martí et al. (marti98 (1998)) using the Very Large Array, but the discovery of radio jets was only possible when the source was observed at milliarcsecond scales with the Very Long Baseline Array (Paredes et al. paredes00 (2000)). The radio emission is persistent, non-thermal and variable, but no strong radio outbursts or periodic variability have been detected (Ribó et al. ribo99 (1999), hereafter R99; Ribó ribo02t (2002)).

In the optical band LS~5039/RX~J1826.2$-$1450 appears as a bright $`V`$=11.2, O6.5V((f)) star showing little variability on timescales of months to years (Clark et al. clark01 (2001)). Variations of $``$0.4 mag have been reported in the infrared ($`H`$ and $`K`$ bands) but no obvious mechanisms for such variability have been proposed (Clark et al. clark01 (2001)). Recent studies show that LS~5039 is a runaway system moving away from the Galactic plane with a total systemic velocity of $``$150 km s-1 and a component perpendicular to the Galactic plane larger than 100 km s-1 (Ribó et al. ribo02 (2002); McSwain & Gies mcswain02 (2002)). The orbit of LS~5039/RX~J1826.2$-$1450 was first studied by McSwain et al. (mcswain01 (2001)), who found that LS~5039/RX~J1826.2$-$1450 is a short-period, $`P_{\mathrm{orb}}=4.117\pm 0.011`$ day and highly eccentric, $`e=0.41\pm 0.05`$, system. Spectroscopic observations carried out in July 2002 by Casares et al. (casares03 (2003)) seem to confirm the orbital period of the system, but neither confirm nor reject the proposed eccentricity.

In the X-rays LS~5039/RX~J1826.2$-$1450 has been studied by R99 using RXTE data. The X-ray timing analysis indicates the absence of pulsed or periodic emission on time scales of 0.02–2000 seconds. The source spectrum is well represented by a power-law model, plus a Gaussian component describing a strong iron line at 6.6 keV. Significant emission is seen up to 30 keV, and no exponential cut-off at high energy is required. The X-ray luminosity has been shown to be compatible with accretion from the stellar wind of the optical companion (McSwain & Gies mcswain02 (2002)).

One of the open questions that remains to be solved is the nature of the compact object of the system as no conclusive data exists. The absence of X-ray pulsations and exponential cut-off at $``$30 keV and the persistent non-thermal emission favours the black-hole classification (R99). On the other hand, the value of the mass function of the system $`f(m)=0.00103\mathrm{M}_{}`$ (McSwain et al. mcswain01 (2001)) and the characteristics of its accretion mechanism (McSwain & Gies mcswain02 (2002)) appear to be more consistent with a neutron star companion.

In this work we present an analysis of the X-ray timing and spectroscopic properties of RX~J1826.2$-$1450 using the narrow field instruments on board the BeppoSAX satellite. This observation permits us to derive a more accurate value of the absorption to the source and study its long-term X-ray variability. We also present optical spectroscopic observations around the H$`\alpha `$ line carried out at the Skinakas Observatory.

## 2 Observations and results

### 2.1 X-ray observations

LS~5039/RX~J1826.2$-$1450 was observed by BeppoSAX on October 8, 2000 (JD 2451825.79–2451826.71) for about 80 ks. BeppoSAX (Boella et al. 1997a ) carried four narrow field instruments (LECS, MECS, HPGSPC and PDS) and a wide field camera (WFC). Here we report on the results from the lower energy detectors: the Low-Energy (LECS) and Medium-Energy (MECS) Concentrator Spectrometers. They consisted of a set of four and three X-ray concentrators respectively together with imaging gas scintillation proportional counter detectors located at the focal planes. In the LECS (Parmar et al. parmar97 (1997)), one of these detectors was sensitive to X-rays in the energy range 0.1–10 keV, while the other three covered an energy range of 1.3–10 keV. The MECS (Boella et al. 1997b ) operated in the energy range 1.6–10 keV. Each detector unit provided a total collecting (geometrical) area of $``$124 cm2, an energy resolution of $`<`$ 8.8 % at 6 keV and a maximum time resolution of 16$`\mu `$s. Note that one of the MECS became inoperative on May 9, 1997 due to a fault in the unit’s gas cell high voltage supply.

Light curves and spectra were extracted from circular regions with a radius of 4′ for both the LECS and the MECS instruments. We selected two source-free regions (of radius 4′ each) in the field of view, and obtained a mean light curve in order to account for the instrumental background. For the spectral analysis, the background was obtained from a blank field using a region similar in size and position to the source extraction region.

The X-ray intensity in the energy range 1.6–10 keV decreased with time from a value of $``$0.06 MECS count s-1 at the beginning of the observation to $``$0.03 MECS count s-1 at the end (Fig. 1). A fit to a constant gave an unacceptable reduced $`\chi ^2`$ ($`>`$2). However, adding a linear term with negative slope we obtained a reduced $`\chi ^2`$1. An F-test gave a probability of less than 0.7% that this improvement in the quality of the fit could occur by chance. The flux variability occurs at higher energies. Below $``$3 keV the flux remains fairly constant (Fig. 1).

Acceptable fits of the X-ray energy spectra were obtained with an absorbed power-law model for both LECS (Fig. 2) and MECS (Fig. 3). Table 1 shows the best-fit spectral parameters. The power law was modified by a low-energy absorption component of the form $`e^{N_\mathrm{H}\sigma (E)}`$, where $`N_\mathrm{H}`$ is the equivalent hydrogen column density and $`\sigma (E)`$ the photoelectric cross-section of Morrison & McCammon (mm83 (1983)). The power-law index and the hydrogen column density have values of $``$1.8 and $``$ 1.0 $`\times `$ 1022 cm-2, respectively (see discussion in Sect. 3.1). The LECS (0.1–5 keV) and MECS (1.6–10 keV) absorbed X-ray luminosity was $`1.9\times 10^{33}`$ erg s-1 and $`4.6\times 10^{33}`$ erg s-1, respectively, assuming a distance of 2.9 kpc (Ribó et al. ribo02 (2002)).

### 2.2 Optical observations

Optical spectroscopic observations were made with the 1.3 m telescope of the Skinakas Observatory in Crete (Greece) on the night of July 26, 1999 and September 10, 2002. The instrumental set-up consisted of a 2000$`\times `$800 ISA SITe CCD and a 1302 l mm-1 grating, giving a dispersion of $``$1 Å/pixel. The reduction of the spectra was made using the STARLINK Figaro package (Shortridge et al. shortridge01 (2001)), while their analysis was performed using the STARLINK Dipso package (Howarth et al. howarth98 (1998)).

The wavelength coverage of the observations included the H$`\alpha `$ line ($`\lambda 6563`$ Å). The strength of the stellar wind emission in H$`\alpha `$ provides information about the mass-loss rate in O-type stars (Puls et al. puls96 (1996) and references therein). The measurement of the equivalent width of the H$`\alpha `$ line (EW(H$`\alpha `$)) in LS~5039 is hampered by the presence of other stellar lines, especially He ii $`\lambda `$6527. Nevertheless, when the emission is low the main factor contributing to the uncertainty in the equivalent width is the difficulty in obtaining the intensity of the continuum. In Table 2 we give the values of EW(H$`\alpha `$) as contemporaneous as possible with the X-ray measurements.

## 3 Discussion

### 3.1 Hydrogen column density

The RXTE observation reported by R99 was not able to provide a reliable value of the absorption to the source due to the insensitivity of the detectors to X-rays below 2.5 keV, where most of the interstellar absorption of X-ray photons takes place. The sensitivity of the BeppoSAX narrow field instruments to X-rays below 2 keV, especially LECS, allows us to derive a more accurate value of the hydrogen column density, $`N_\mathrm{H}=1.0_{0.3}^{+0.4}\times 10^{22}`$ cm-2. Assuming $`E(BV)=A_V/3.1`$ (Rieke & Lebofsky rieke85 (1985)) the X-ray absorption can be converted into optical extinction by means of the empirical relations of Predehl & Schmitt (predehl95 (1995)), to obtain $`E(BV)=1.8_{0.5}^{+0.7}`$, or those of Ryter et al. (ryter75 (1975)) and Gorenstein (gorenstein75 (1975)), to obtain $`E(BV)=1.5_{0.4}^{+0.6}`$. These values of $`E(BV)`$ are considerably higher than the unreliable value of $`E(BV)0.4`$ derived from the PCA/RXTE observation (R99).

On the other hand, optical photometry of LS~5039 can be found in Drilling (drilling91 (1991)) and Lahulla & Hilton (lahulla92 (1992)) who gave a value for the colour $`(BV)=0.95\pm 0.01`$ and Clark et al. (clark01 (2001)), who reported $`(BV)=0.85\pm 0.02`$. Using the average of these values and assuming an intrinsic colour index of $`(BV)_0=0.30\pm 0.02`$ for an O6.5V star (Lejeune & Scharer lejeune01 (2001)) we obtain a colour excess of $`E(BV)=1.2\pm 0.1`$. The extinction obtained from the BeppoSAX observation is somehow higher than (although statistically consistent with) that derived from the optical observations. One explanation for this difference might be the fact that the hydrogen column density is more sensitive to absorption by both interstellar and intrinsic (i.e. in the vicinity of the X-ray source) matter, while the optical indicators would mainly account for the interstellar absorption. The presence of significant amounts of cold matter in the surroundings of the system is supported by the detection of the relatively strong iron (equivalent width of 0.75$`\pm `$0.06 keV) emission line at 6.6 keV in the 3–30 keV X-ray spectrum of LS~5039/RX~J1826.2$-$1450 (R99). Unfortunately, due to limited signal-to-noise ratio no line is detected in the BeppoSAX data.

### 3.2 Absence of X-ray eclipse

Using the ephemeris of McSwain et al. (mcswain01 (2001)), the BeppoSAX observation covers the orbital phase from 0.89 to 0.11. Note that this observation took place just four days after the periastron date used by these authors in their ephemeris, i.e., JD 2,451,822.12$`\pm `$0.09. Consequently, our determination of the orbital phase is not strongly affected by the lack of accuracy in the value of the orbital period ($`P_{\mathrm{orb}}=4.117\pm 0.011`$ days). Hence, periastron took place on JD 2,451,826.24$`\pm `$0.09, or at $`4.2\pm 0.8\times 10^4`$ s in Fig. 1. The epoch of inferior conjunction of the optical star, which would correspond to the time of an X-ray eclipse, is 0.25 days after periastron or at orbital phase 0.06 (McSwain et al. mcswain01 (2001)), which corresponds to JD 2,451,826.49$`\pm `$0.09, or to $`6.4\pm 0.8\times 10^4`$ s in Fig. 1. The duration of a possible eclipse would last less than 0.2 days ($`1.7\times 10^4`$ s) for $`M_{\mathrm{opt}}=40\mathrm{M}_{}`$, $`R_{\mathrm{opt}}=10\mathrm{R}_{}`$ and $`i=90\mathrm{°}`$ (shorter for a lower inclination and an eventual lower mass for the primary) and hence it should have been detected by BeppoSAX.

Although the light curve does not show statistical significant evidence for an eclipse, the variability increases between periastron and inferior conjunction, possibly indicating some kind of distortion of the stellar wind when the compact star finds itself at minimum distance from the primary. There is also a trend for the count rate to decrease as the source approaches inferior conjunction. This trend is detected at energies above 3 keV, but not at lower energies, contrary to what would be expected if this was due to a change in absorption, as it happens in Cygnus~X-1 (Wen et al. wen99 (1999)). Therefore, this decrease in the hard X-ray flux during the BeppoSAX observation is probably intrinsic, and may reflect changes in the accretion rate due to orbital motion.

The lack of eclipses allows us to set an upper limit on the inclination of the system $`i`$. This implies, in the case of a spherical companion (no photometric ellipsoidal modulation has been found to date), that $`\mathrm{cos}i>R_{\mathrm{opt}}/r`$, where $`R_{\mathrm{opt}}`$ is the radius of the optical companion an $`r`$ is its distance to the compact object at inferior conjunction (phase 0.06). Assuming a mass of 40$`\mathrm{M}_{}`$ and a radius of 10$`\mathrm{R}_{}`$ for an O6.5V star (Howarth & Prinja howarth89 (1989)), the canonical mass for a neutron star (1.4$`\mathrm{M}_{}`$) and the orbital parameters of McSwain et al. (mcswain01 (2001)), $`e=0.41\pm 0.05`$ and $`P_{\mathrm{orb}}=4.117\pm 0.011`$ days, we obtain an upper limit for the inclination of the system of $`i<66\pm 2\mathrm{°}`$, where the upper and lower values correspond to the range of eccentricities given above. This inclination is consistent with the one, $`i_f`$, derived from the mass function, $`f(m)=0.00103\pm 0.00020\mathrm{M}_{}`$ (McSwain et al. mcswain01 (2001)), for the same parameters, namely, $`i_f=60\mathrm{°}`$.

Note that OB stars in X-ray binary systems may be less massive than the corresponding field stars of the same spectral type (Kaper kaper01 (2001)). A decrease of the mass of the optical companion by about 10% reduces $`i_f`$ to $`54\mathrm{°}`$, while the upper limit hardly changes. On the other hand, assuming $`M_{\mathrm{opt}}=40\mathrm{M}_{}`$ and the above mentioned orbital parameters Ribó et al. (ribo02 (2002)) derived an upper limit for the mass of the compact object of 9 $`\mathrm{M}_{}`$ (in order to avoid breackup speed of the primary once $`V_{\mathrm{rot}}\mathrm{sin}i`$ is known). For a $`9\mathrm{M}_{}`$ black-hole companion a similar analysis yields $`i<68\pm 2\mathrm{°}`$ and $`i_f=8.5\mathrm{°}`$. Similar values are obtained if we decrease the mass of the optical companion by about 10%. In either case, the lack of X-ray eclipses constrains the inclination of the system to values smaller than 70°, which is clearly compatible with those found using the mass function of the system.

### 3.3 Long-term X-ray variability

In order to study the long-term X-ray variability of LS~5039/RX~J1826.2$-$1450 we obtained the X-ray luminosity in the same energy ranges as given in previous missions. The value of the ROSAT 1996 October 0.1–2.4 keV (Motch et al. motch97 (1997)) and the RXTE 1998 February 3–30 keV (R99) X-ray luminosities are 7.1 $`\times `$ 1033 erg s-1 and 5.2 $`\times `$ 1034 erg s-1, respectively. The 0.1–2.4 keV X-ray luminosity obtained with LECS is 6.6 $`\times `$ 1032 erg s-1 and that obtained with MECS in the energy range 3–30 keV is 8.0 $`\times `$ 1033 erg s-1 (derived assuming the model given in Table 1). This observation took place in 2000 October. All values are for an assumed distance of 2.9 kpc (Ribó et al. ribo02 (2002)). Clearly the source was in a low-flux state during the BeppoSAX observation (about one order of magnitude fainter) when compared to the epoch of the ROSAT and RXTE observations. As explained above this low flux, however, cannot be accounted for by invoking a partial X-ray eclipse. We have also inspected archival X-ray data and found an ASCA observation of LS~5039/RX~J1826.2$-$1450 made on 1999 October 4. We have obtained the X-ray flux in the energy range 0.3–10 keV, shown in Table 2, which is between the RXTE and BeppoSAX measurements. In addition, a similar X-ray flux has been recently detected with CHANDRA (Miller, private communication), and is also reported in Table 2.

We also studied the long-term optical variability of the primary by measuring the strength of the H$`\alpha `$ line. Table 2 shows the values of the H$`\alpha `$ equivalent width together with the X-ray flux from the missions that have observed LS~5039/RX~J1826.2$-$1450. The H$`\alpha `$ emission was weakest during the BeppoSAX observation in October 2000 and strongest a few months after the RXTE pointing in February 1998, when the X-ray emission was also strong (since the line is in absorption, the smaller $`EW`$(H$`\alpha `$) the stronger the emission). Moreover, very similar $`EW`$(H$`\alpha `$) are found contemporaneous to the ASCA and CHANDRA observations, as the detected X-ray fluxes are, and also between the $`EW`$(H$`\alpha `$) measurements contemporaneous to the RXTE and BeppoSAX observations. Hence, there is a trend in the sense that the lower the emission in H$`\alpha `$, the lower the observed X-ray flux. All these results suggest that a decrease in the mass-loss rate from the primary induces a decrease in the emitted X-ray flux of the binary system. In this context, we propose that changes in the mass-loss rate from the primary, reflected in the change of $`EW`$(H$`\alpha `$), can account for a change in the mass accretion rate onto the compact object, and hence explain the detected changes in X-ray luminosity.

On the other hand, since this binary system has a high eccentricity, we cannot rule out a priori that the detected variations in the X-ray flux are not due to orbital variability. If the mass transfer in LS~5039/RX~J1826.2$-$1450 occurs via the strong stellar wind of the optical primary then, given the high eccentricity of the system, some variations in the X-ray luminosity with the orbital phase should be expected. The fraction of the stellar wind captured by the compact object depends on its velocity relative to the wind. Also, variations in the density of the wind may produce changes in the amount of absorption material throughout the orbit giving rise to the X-ray flux variations. Orbital flux variations have been seen in many high-mass X-ray binaries powered by stellar wind accretion: 4U~1700$-$37 (Haberl & Day haberl92 (1992)), GX~301$-$2 (Haberl haberl91 (1991)), Vela~X-1 (Haberl & White haberl90 (1990)), 4U~1907+09 (Marshall & Ricketts marshall80 (1980)). With its small orbit and early-type optical primary LS~5039/RX~J1826.2$-$1450 is expected to share many similarities with these systems.

However, a purely Bondi-Hoyle accretion model (Bondi bondi52 (1952)) does not reproduce the observations. Figure 4 shows the expected variability of the X-ray luminosity in LS~5039/RX~J1826.2$-$1450 using a spherically symmetrical Bondi-Hoyle accretion model. For each orbital phase and using the formulae given in the Appendix, we computed the binary separation, the orbital and wind velocities, the relative velocity, the wind density, the accretion rate and finally the X-ray luminosity. We note that this is a bolometric unabsorbed X-ray luminosity, to be compared to the much lower observed luminosity in the range 1.6–10 keV (top panel of Fig. 1). Although this discrepancy could also be caused by the lack of precise knowledge of the parameter values used in the model, there are two strong arguments against this model. First, we should expect a variation in flux of at least a factor of 2 within the phase interval of the BeppoSAX observations as follows: starting with a gradual rise, peaking smoothly around phase 0.02 ($`5\times 10^4`$ s in Fig. 1) and ending with a gradual decrease. Second, since the BeppoSAX observation took place during periastron the X-ray luminosity should be the largest, contrary to what it is observed (see Table 2).

Although it is clear that detailed observations covering a full orbital period are needed in order to better study and model the orbital X-ray variability of LS~5039, it seems unlikely that it can be responsible for the variations of one order of magnitude quoted in Table 2. Therefore, we conclude that this variability may be due to changes in the stellar wind mass-loss rate on longer timescales.

## 4 Conclusions

The main conclusions of the analysis of our BeppoSAX observation of the microquasar LS~5039/RX~J1826.2$-$1450 can be summarized as follows:

1. We have obtained the most accurate value of the hydrogen column density to date: $`N_\mathrm{H}=1.0_{0.3}^{+0.4}\times 10^{22}`$ cm-2. The derived colour excess is compatible but higher than the one obtained from optical photometry alone, suggesting that there might be intrinsic absorption within the binary system orbit.
2. Our observation took place near inferior conjunction, and no X-ray eclipse has been detected. This fact has allowed us to compute upper limits for the inclination of the binary system for different masses of the compact object: $`i<66\mathrm{°}`$ for $`M_\mathrm{X}=1.4\mathrm{M}_{}`$ and $`i<68\mathrm{°}`$ for $`M_\mathrm{X}=9\mathrm{M}_{}`$. These values are compatible with the ones obtained using the mass function.
3. We have found hints of a correlation between the measured $`EW`$(H$`\alpha `$) and the observed X-ray flux on timescales of years, in the sense that the stronger the optical emission the higher the X-ray flux. In this context, we have suggested that a decrease in the mass-loss rate from the primary induces a decrease in the mass accretion rate onto the compact object, and hence a decrease in the observed X-ray flux.
4. A simple Bondi-Hoyle accretion model (without taking into account a varying absorption) cannot explain the X-ray variability because it predicts an increase in the X-ray flux near periastron that is not observed. This may happen because the model is very sensitive to some parameters not yet very well constrained, like the eccentricity and the wind velocity. However, the fact that no significant variability is detected near periastron suggests that the observed variability in the X-ray flux on timescales of years is due to changes in the mass-loss rate of the primary. An X-ray monitoring during a full $``$4 d orbital cycle would clarify these ideas.

## Appendix: Orbital Bondi-Hoyle X-ray variability

An order of magnitude estimate of the orbital X-ray variability can be obtained using the Bondi-Hoyle accretion model. The accretion luminosity is then computed by means of the simple formula:

$$L_\mathrm{X}\frac{GM_\mathrm{X}\dot{M}}{R_\mathrm{X}},$$
(1)
where $`G`$ is the gravitational constant, $`M_\mathrm{X}`$ and $`R_\mathrm{X}`$ the mass and radius of the compact object and $`\dot{M}`$ its variable accretion rate. Assuming a neutron star, we will proceed adopting $`M_\mathrm{X}=1.4\mathrm{M}_{}`$ and $`R_\mathrm{X}=10`$ km. If the wind accretion onto the compact companion can be described as a spherical inflow (Bondi bondi52 (1952)), the accretion rate as a function of the binary separation $`r`$ is given by:

$$\dot{M}=4\pi (GM_\mathrm{X})^2\frac{v_\mathrm{w}(r)\rho (r)}{v_{\mathrm{rel}}^4(r)},$$
(2)
where $`v_\mathrm{w}(r)`$ is the wind velocity, $`\rho (r)`$ is the density of the stellar wind and $`v_{\mathrm{rel}}(r)`$ its supersonical velocity relative to the neutron star. We will approximate the wind velocity $`v_\mathrm{w}(r)`$ by a typical $`\beta `$-law appropriate for O-type stars:

$$v_\mathrm{w}(r)=v_{\mathrm{}}\left[1\frac{R_{\mathrm{opt}}}{r}\right]^\beta .$$
(3)
Here, $`v_{\mathrm{}}`$ is the wind velocity at infinity and $`\beta `$ a power index. Assuming a simple radial flow, the corresponding wind density can be derived from the continuity equation as:

$$\rho (r)=\frac{\dot{M}_{\mathrm{opt}}}{4\pi r^2v_\mathrm{w}(r)}.$$
(4)
Following McSwain & Gies (mcswain02 (2002)), suitable values for LS~5039 are $`v_{\mathrm{}}=2850`$ km s-1, $`\beta =0.8`$ and $`\dot{M}_{\mathrm{opt}}=8\times 10^7\mathrm{M}_{}`$ yr-1. To proceed, we have solved numerically the Keplerian motion of the neutron star around the $`M_{\mathrm{opt}}=40\mathrm{M}_{}`$ companion calculating both the eccentric and true anomalies ($`E`$ and $`V`$). The relative velocity between the wind and the neutron star is then given by:

$`v_{\mathrm{rel}}^2(r)=\left[{\displaystyle \frac{na^2}{r}}\mathrm{sin}E+v_\mathrm{w}(r)\mathrm{cos}V\right]^2+`$ (5)
$`+\left[{\displaystyle \frac{na^2\sqrt{1e^2}}{r}}\mathrm{cos}Ev_\mathrm{w}(r)\mathrm{sin}V\right]^2,`$
where $`a`$ is the semimajor axis and $`n=2\pi /P_{\mathrm{orb}}`$ is the mean angular motion.

###### Acknowledgements.

We are grateful to J. Miller for permission to include in this paper unpublished results from his CHANDRA observation. We thank useful comments and suggestions from an anonymous referee. P. R. is a researcher of the programme Ramón y Cajal funded by the University of Valencia and Spanish Ministery of Science and Technology. M. R., J. M. P. and J. M. acknowledge partial support by DGI of the Ministerio de Ciencia y Tecnología (Spain) under grant AYA2001-3092, as well as partial support by the European Regional Development Fund (ERDF/FEDER). During this work, M. R. has been supported by a fellowship from CIRIT (Generalitat de Catalunya, ref. 1999 FI 00199). J. M. has been aided in this work by an Henri Chrétien International Research Grant administered by the American Astronomical Society, and has been partially supported by the Plan Andaluz de Investigación of the Junta de Andalucía (ref. FQM322). Skinakas Observatory is a collaborative project of the University of Crete, the Foundation for Research and Technology-Hellas and the Max-Planck-Institut für Extraterrestrische Physik.


text:
# The envelope of IRC+10216 reflecting the galactic light Based on observations made with the Antu 8-m VLT (ESO; program 63I-0177A), with the Canada France Hawaii 3.6-m Telescope (CNRS, NRC, UH) and the 1.20-m telescope of Haute-Provence Observatory (CNRS)

## 1 Introduction

Stars of intermediate initial mass can have very high mass-loss rates when evolving on the asymptotic giant branch (AGB). This mass loss produces extended, dusty, optically thick envelopes that strongly mask the central star in the visible domain, possibly leaving no optical counterpart. However, these envelopes are externally illuminated by the stars of the galactic disk. If deep enough, a direct image of the envelope can show the reflection of this external photon source by circumstellar grains. Two such cases of AGB reflection nebulae are presently known. The best case is the nearby carbon star IRC+10216, with an exceptionally massive circumstellar envelope that is clearly seen in the ambient galactic light in the $`B`$ and $`V`$ bands (Mauron & Huggins mauron99 (1999)). The envelope is detected up to about 3 from center, and appears fairly round, which is qualitatively consistent with a spherically symmetric dust shell and an isotropic galactic field (at least in $`B`$ and $`V`$). There is also a central plateau of diameter $``$ 20′′, (i.e. no central peak), which is prominent in $`B`$ and proves that the source of photons is not the central star. The second case is CRL 3116, another carbon AGB object at about 1 kpc from the Sun (Crabtree & Rogers crabtree92 (1992)), but this object is much less studied than IRC+10216.

In the present paper, we follow the works of Martin & Rogers (martin87 (1987), hereafter MR87) and Crabtree & Rogers (crabtree92 (1992)), and our goal is to make a quantitative investigation of these reflection nebulae in order to derive information on the interstellar radiation field (ISRF) and/or the circumstellar dust properties. Briefly, the ISRF is due to the integrated light of galactic disk stars and a diffuse galactic light component. It is an important quantity involved in many problems concerning the interstellar medium or circumstellar shells, such as grain heating, their far-infrared emission and photochemistry (e.g., Mathis, Mezger & Panagia mathis83 (1983), van Dishoek vandishoek94 (1994), Glassgold glassgold96 (1996)). The intensity of the local ISRF (measured from Earth) is relatively well known in the visible range, but the situation is more controversial in the ultraviolet, with differences between authors of at least a factor of 2 at 1500Å (see Fig. 2 of van Dishoek vandishoek94 (1994)). Concerning the spatial ISRF variations within the solar neighbourhood, one expects rather small variations, except for peaks close to very luminous stars, especially in the far ultraviolet region (Habing habing68 (1968), Jura jura74 (1974)).

Martin & Rogers (martin87 (1987)) found that in order to account for the $`V`$ surface brightness of the IRC+10216 plateau, which was relatively faint compared to their model, a possible solution was to assume an ISRF smaller than that of Mathis et al. (mathis83 (1983)) by a factor of $``$ 2.5. They suggested that a faint ISRF could be due to the height of this object above the galactic plane, which is $`z=90`$ pc for an assumed distance $`d=130`$ pc. However, according to Mattila (1980a ), the dependance on $`z`$ is predicted to be much smaller than a factor of 2, at least in the optical domain. Another indication that the ISRF may have no strong minima comes from the COBE observations of interstellar dust at high galactic latitude: the narrow range of dust temperature suggests ISRF variations smaller than 30% (Lagache et al. lagache98 (1998)).

Martin & Rogers also noted two other possible explanations for the low plateau $`V`$-band brightness: the first is to assume a grain albedo lower than originally adopted in their model; the second is to envisage an even smaller grain radius than assumed in their model (0.05 $`\mu `$m) and to increase the mass loss rate in the external layers. However more recent studies of the dust of IRC+10216 point to significantly larger grain sizes: a dominant radius of 0.16 $`\mu `$m is found by Groenewegen (groenewegen97 (1997), GR97 hereafter). Distribution sizes including grains as large as 0.2 or 0.4 $`\mu `$m are also favoured by Ivezic & Elitzur (ivezic96 (1996)) and by Skinner, Meixner & Bobrowsky (skinner98 (1998)).

These considerations led us to obtain additional images of the ISRF-reflecting nebula around IRC+10216. Our goal is to reanalyse the intensity, colors, size, and shape of this nebula. In particular, this paper presents new surface brightness measurements in the three $`UBV`$ bandwidths (a photometric sequence of field stars is given in appendix A). Details about the observations and the data reduction are given in Sect. 2, and the results are presented in Sect. 3. A specific Monte-Carlo radiative transfer code was developped for this program and is described in Sect. 4, with complementary considerations in Appendix B. The data are quantitatively interpreted and discussed in Sect. 5, and a summary of this work is presented in the conclusion.

## 2 Observations and reductions

The first set of observations was carried out at Haute-Provence Observatory with the 1.20-m telescope, equipped at the Newtonian focus with a CCD camera for $`UBVRI`$ imaging. The chip is a thin back-illuminated 10242 Tektronix, with a read-out noise of 8.5 $`e^{}`$ and a gain of 3.5 $`e^{}`$ ADU-1. This instrument provides a field of 11.7$`\times `$11.7 with a scale of 0.6845 arcsec pixel-1. More details can be found in Ilovaisky (ilovaisky97 (1997)). The nights were clear but the seeing was 2′′. Flat-fields were done on the dome and on the sky. The $`UBV`$ images of IRC+10216 obtained are listed in Table 1. Several images were also taken each night on the calibration fields NGC 2403 (Zickgraf et al. zickgraf90 (1990)) and M13 (Arp & Johnson arp55 (1955)). All CCD frames were reduced in the standard way, with bias substraction and division by averaged flat-fields. These data permitted the establishment of a $`UBV`$ photometric sequence containing 16 stars in the field of IRC+10216 (details in Appendix A).

The CFHT observations are listed in Table 1 and were described by Mauron & Huggins (mauron99 (1999)). They resulted into two deep $`B`$ and $`V`$ images, with a field of 8.7$`\times `$8.7 and a scale of 0.436 arcsec pixel-1. Because no photometric standards were observed at CFHT, calibration was achieved by considering unsaturated stars of the OHP sequence.

The VLT observations were done in service mode with the Antu 8-m telescope, equipped with the FORS1 focal reducer (Table 1). The detector is a 20482 thinned 24 $`\mu `$m pixel Tektronix chip. The field of view is 6.83$`\times `$6.83 with a pixel size of 0.200 arcsec pixel-1. More details on the instrumentation can be found in the FORS1 Manual at ESO Web site. The usable data consists of three 8-min $`U`$ frames and two 2-min $`V`$ frames. A number of other frames in $`UBV`$ were acquired on different nights, but had to be rejected for this study due to an excessively bright sky or poor photometric conditions. The frames were merged into one $`U`$ and one $`V`$ flat-fielded image \[This 4 mn VLT $`V`$-band image is much less deep than the CFHT one presented in Mauron & Huggins (mauron99 (1999)) and is not shown here; for the $`U`$-band image, see below\]. During the nights of observation, one or two exposures per filter were supplied of calibration fields (Mark A, SA 101 and PG 0942-029, detailed on the FORS ESO Web site). These exposures did not allow an independant determination of the atmospheric extinction coefficients of these nights, but if plausible values are assumed ($`k_U`$=0.62, $`k_V`$=0.14 mag/airmass), one finds that the $`U`$ and $`V`$ photometry of the field stars agree very well with the OHP results, with differences of 0.04 mag in $`V`$, and 0.03 mag in $`UV`$. Therefore, all measurements on the nebula were calibrated using the OHP $`UBV`$ sequence.

After bias subtraction and flat-field correction, our goal was to obtain surface photometry of the plateau in the $`UBV`$ bands, and to derive the surface brightness profiles. The main difficulty is to estimate the background sky level to be subtracted, for several reasons. First, the plateau is intrinsically faint, about 4% of the sky brightness in $`V`$. This implies that, if an accuracy of $``$ 10% is desired for the plateau photometry, the sky brightness far from the nebula has to be determined with a relative error of $``$ 0.4%. The flat field correction should also be achieved with this accuracy over large or medium angular scales, but in practice the corrected images show small but non-negligible gradients. These gradients can be corrected in part by dividing by second order polynomial functions, but residual variations remain at intermediate scales. The situation is also complicated by the fact that there are three bright stars in the field: they produce haloes, various ghosts, and significant diffuse parasite light that is not uniform. Nevertheless, these defects depend on the intrumentation used, and we could perform cross-verifications between OHP, CFHT, and VLT data. The surface photometry of the nebula was obtained with azimuthally averaged radial profiles. In the $`U`$-band, the center position given by the $`V`$ frames was adopted . The profiles are achieved after masking all stars, galaxies, bright stars and their haloes, and other obvious defects (see MH99 for more details).

## 3 Observational results

The reduction of all images resulted in the final plateau $`UBV`$ surface photometry given in Table 2. The $`UBV`$ measurements refer to the top of the plateau, i.e. at $`r=10`$′′, on azimutally averaged profiles. The presence of a central nearly point-like peak due to radiation from the star in the CFHT and VLT $`V`$ images was taken into account and a small correction was done when estimating the plateau maximum intensity at $`r=10`$′′. No $`V`$ plateau photometry could be achieved from the OHP data because of an excessively strong central peak. The $`B`$ and $`V`$ surface brightnesses of Table 2 are not very different to previous values given by MH99, and supersede them. In $`V`$, the plateau is found to be 0.24 mag brighter than that found by Crabtree et al. (crabtree87 (1987)).

Figure 1 displays a composite image of the envelope in the $`U`$ band, made by summing all the long OHP and short VLT $`U`$ exposures (the summed number of photons per arcsec2 emitted by the nebula is comparable in OHP and VLT observations). Comparison with the $`B`$ and $`V`$ images displayed in MH99 shows that, in the $`U`$ filter, the nebula may be slightly brighter on the South-west side, which is the direction towards the galactic equator. This might be due to the UV-bright OB stars concentrated in the galactic plane, but a deeper image is obviously necessary to confirm this asymmetry.

Figure 2 shows azimuthally averaged profiles of the nebula. Each point corresponds to the signal average over a 4′′-wide annulus. No sky subtraction has been done here. The ordinates are CCD counts (ADU) and the scales have been adjusted so that both the plateau and the sky are located at about the same level (dashed line) for any instrument and filter. This scaling could not be done for OHP $`V`$ because the plateau is not well detected, and this profile is displayed here only to show the nebula extension at faint levels. This figure shows that there is no noticeable difference in the nebula shape as a function of wavelength, given the experimental uncertainties. It is worth noticing that despite the very low net signals in the $`U`$-band, i.e. 1.5 ADU at VLT, 5 ADU at OHP, the $`U`$ profiles are well recovered thanks to the averaging over a large number (thousands) of pixels located in the annuli. The profiles show that the nebula is detected up to $``$ 100 or $``$ 200 arcsec, depending on instrument and bandwidth. These differences are instrumental effects. They reflect both various depths of the exposures and lack of background flatness at very faint levels.

The depth of the central depression is best measured in the $`B`$ band. In this band the central core of IRC+10216 is a very faint point-like source, with $`B25.2`$ in the CFH frame. Its contribution can be easily subtracted, and the depression depth is found to be 17$`\pm `$3 percent in $`B`$. Estimating the depth in $`U`$ is much less accurate due to poor S/N in the central bin ($`r<4`$′′).

The half maximum radius of the profile is also best measured in $`B`$ and is $`R_{1/2}=29`$′′ $`\pm `$ 2′′ from OHP and CFH data. There is no evidence for a larger $`R_{1/2}`$ value in $`U`$, taking into account experimental uncertainties (Fig. 2). Finally in the $`V`$ band, $`R_{1/2}`$ is measurable with less accuracy, and is estimated to be around 27′′.

## 4 Modelling

In order to interpret the observational data, we have to adopt quantitative characteristics for the ISRF and the circumstellar dust envelope (CSE), and the scattering of light by the dust has to be calculated. Our approach to the ISRF and the CSE parameters is to consider the most recent and plausible values derived from other studies (although as seen below, significant uncertainties exist), and to investigate whether they can explain the observations.

### 4.1 Parameters for the ISRF

We have considered two different sources for the ISRF. The first is that of Mathis, Metzger & Panagia (mathis83 (1983), MMP83), who carried out a synthesis of various measurements and models in order to derive ISRF values as a function of galactocentric distance, from far UV to far IR wavelengths. Table 3 reproduces the $`\lambda `$, $`f(\lambda )`$ values of their Table A3 for the solar neighbourhood (i.e. for $`D_G=10`$ kpc as adopted by them). Corresponding apparent $`UBV`$ magnitudes and colors were derived with the appropriate formulae of Allen (allen73 (1973)), which are indicated in the Table 3 notes.

The second work we considered is that of Mattila (1980a ) who calculated synthetic spectra of the integrated stellar light (ISL) between 3 000 Å and 10 000 Å in the solar neighbourhood. The diffuse galactic light, another component of the ISRF in addition to direct starlight and which comes from light scattered by dust, is not modelled in that work. However, according to Mattila (1980b ), its effect on the total radiation density averaged over the whole sky is estimated to be small, of the order of 20%. Here, we have chosen to ignore this small component for clarity. Table V of Mattila (1980a ) provides $`UBV`$ surface brightness estimates of the ISL, assumed here to be equal to the ISRF, in the solar neighbourhood and in units of S10, where one S10 is one 10th magnitude star per square degree. They were converted into apparent $`UBV`$ magnitudes by integrating over the sky (41253 sq. deg.). The corresponding fluxes $`f(\lambda )`$ derived with Allen’s formulae are also given in Table 3.

Table 3 shows that these two ISRFs agree in $`B`$ and $`V`$ ISRF fluxes (and $`B`$-$`V`$ color), but a difference of 0.45 mag. is found in the $`U`$ band (i.e. a factor of 1.5), in the sense that the Mathis ISRF is bluer in $`U`$-$`B`$. Mattila (1980a , 1980b ) found in his synthetized ISL a large flux jump between 3 500 Å and 4 200 Å which results from the Balmer jump of A-type stars and the discontinuity at 4 000 Å from late-type stars. In contrast, the MMP83 ISRF is expressed by 3 blackbodies $`plus`$ an UV component for $`\lambda <`$ 2 500 Å, and this approximation is said to be within 15 percent (or better) of the various observational data over the full spectral range (see MMP83, and also Mezger et al. mezger82 (1982), their Appendix C). Therefore, at this stage, we find it interesting to keep both approaches for further discussion.

### 4.2 Parameters for the dust envelope of IRC+10216

Concerning the CSE of IRC+10216, we have adopted as a baseline the results of Groenewegen (groenewegen97 (1997), herafter G97), who considered a large number of constraining observations (that probe mostly, but not exclusively, the inner layers of the CSE). The adopted parameters are listed in Table 4. Although the CSE is known to show multiple shells (MH99, Mauron & Huggins mauron00 (2000)), we assume a smooth and spherically symmetric density distribution and will attempt to reproduce azimuthally-averaged profiles, leaving the consequences of the non-radial envelope structure for future study. With an assumed distance of 135 pc, G97 estimates a dust mass-loss rate of $`\mathrm{1.1\; 10}^7`$ M yr$`^1`$beyond 4.47 stellar radii (inner envelope radius). Furthermore, we considered an outer envelope radius of 210” ($`R_{\mathrm{out}}`$ $`=`$ 6 000 $`R_{}`$) corresponding to the most distant material detected by Mauron & Huggins (mauron99 (1999)). This is certainly a lower limit for the IRC+10216 envelope outer radius since Young et al. (young93 (1993)) found dust twice as far away using IRAS data. However, we checked that our results concerning the model brightnesses and colors of the IRC+10216 nebula are independent of $`R_{\mathrm{out}}`$. Considering a larger outer radius would only lead to a slightly larger width of the brightness profiles (see the discussion below).

As recommended by G97, the grains considered here are spherical particles with a single size ($`a`$ = 0.16 $`\mu m`$), and their composition is that of amorphous carbon with optical properties given by Rouleau & Martin (rouleau91 (1991), labelled AC1 in their Table 1; see also Appendix C). This table also indicates another possible choice for the grain radius, $`a`$ = 0.05 $`\mu m`$ (as found by MR87 and discussed below). It has to be noted that silicon carbide grains are also found in the envelope of IRC+10216. However, since (i) only a few percent of the dust grains are found in SiC form, and (ii) their opacity effects are mostly seen around 11 $`\mu `$m and are completely negligible in the optical, we do not consider such grains in our model.

With the parameters of Table 4 (standard model, with $`a=0.16\mu m`$), the radial optical thickness of the envelope in $`B`$ (taken from the inner to the outer radius of the envelope) $`\tau _{B,\mathrm{rad}}`$ is equal to 84.3, in excellent agreement with G97 if he neglected SiC grains. The observed spectral energy distribution of IRC+10216 is well reproduced by our model except for some small departures around 11 $`\mu `$m in the SiC signatures.

The optical thickness of the nebula along the line of sight at $`\theta `$ arcseconds from the center can be written as:

$$\tau _\lambda (\theta )=\pi \tau _{\lambda ,\mathrm{rad}}\frac{\theta _{\mathrm{in}}}{\theta }q[\frac{2}{\pi }arctan(\frac{\sqrt{1u^2}}{u})].$$
(1)
Here, $`\tau _{\lambda ,\mathrm{rad}}`$ is the radial optical thickness of the envelope (i.e. from the center of the nebula up to $`R_{\mathrm{out}}`$); $`\theta _{\mathrm{in}}=R_{\mathrm{in}}/d`$ corresponds to the angular inner radius of the envelope; $`\theta _{\mathrm{out}}=R_{\mathrm{out}}/d`$. The factor $`q`$ is very close to 1.0 and can be ignored: rigourously $`q=1/(1R_{\mathrm{in}}/R_{\mathrm{out}})`$. The quantity $`u`$ is $`\theta /\theta _{\mathrm{out}}`$, and the factor $`[\frac{2}{\pi }arctan()]`$ is the effect of the outer envelope limit, and is equal to 1 if no outer limit is assumed.

With the parameters of Table 4, and for an angular distance from the center of 50′′, the standard model predicts in the $`B`$ filter a line of sight opacity $`\tau _B(50^{\prime \prime })=0.70`$. The factor in brackets due to the external edge is 0.85 and an unlimited envelope with the same $`\tau _{\lambda ,\mathrm{rad}}`$ would have an opacity 15 percent larger, i.e. 0.82 . At 100′′, one finds $`\tau _B(100^{\prime \prime })=0.28`$. Therefore, the model envelope is optically thick in the observed region for offsets between 0 and $``$ 100-200 ′′, and a radiative transfer code is needed to predict the nebula optical surface brightness.

### 4.3 The Monte-Carlo radiative transfer code

In order to model ISRF illumination, we have developed a specific radiative transfer code based on an upgraded version of the one described in Lopez et al. (lopez95 (1995)). Although the code has been built to allow for non-spherical dust density distributions, we consider in the present work spherically symmetric envelopes with a radial power-law distribution of the form $`N(r)=N_{in}\times (R_{in}/r)^n`$, where $`r`$ is the radial distance from the star with $`R_{in}<r<R_{out}`$, and $`N(r)`$ is the number density of grains. The dust grains are assumed to be homogeneous and spherical with a single radius $`a`$. Their optical properties, i.e. the extinction $`Q_{\mathrm{ext}}`$, scattering $`Q_{\mathrm{scat}}`$, and absorption $`Q_{\mathrm{abs}}`$ efficiencies, and the scattering phase function, are derived from optical indices using Mie theory. Finally, the optical thickness of the envelope at the wavelength $`\lambda `$ is defined by

$$\tau _{\lambda ,\mathrm{rad}}=_{R_{in}}^{R_{out}}N(r)\pi a^2Q_{\mathrm{ext}}(\lambda )𝑑r$$
(2)
Two sources of illuminating radiation can be considered in the code: the first is the central star with radius $`R_{}`$ and assumed to radiate as a blackbody of a given effective temperature. This source is considered only for computing the spectral energy distribution (SED), found to be in full agreement with that of G97. The transfer of radiation in the dust and the radiative equilibrium of the grains are numerically solved by a Monte-Carlo method at 30 different wavelengths, and more details are given in Lopez et al. (lopez95 (1995)).

The second source is the ISRF, whose spectral energy distribution is taken either from MMP83 or from Mattila (1980a ). The ISRF is synthesized using a large light-emitting sphere centered on the nebula. The surface brightness of the sphere varies as cos$`\varphi `$, where $`\varphi `$ is the inclination angle to its normal, and is equal for $`\varphi =0`$ to $`f(\lambda )/4\pi `$ (in erg cm2 s-1 Å-1 sr-1), where $`f(\lambda )`$ is the ISRF value as given in Table 3. The radius of this ISRF sphere is larger than the outer radius of the envelope, and it can be shown that each point inside this ISRF sphere receives in the absence of matter a homogeneous and isotropic photon flux equal to the ISRF.

In order to obtain a model image to compare with the $`UBV`$ observations described above, we consider only ISRF photons emitted by the sphere that have been actually scattered by dust grains in the envelope, and we ignore those directly emitted towards the observer and passing through the envelope without any interaction. This has to be done in order to avoid a bright background, which is not observed. The reason is that in the reality (but not in our model), the ISRF photons come mainly from individual point-like sources (the stars). The contribution of diffuse light scattered by the galactic cirrus and interstellar clouds is small, and this is especially true for ISRF photons coming from behind IRC+10216, which is at $`l=221\mathrm{°}`$, $`b=+45\mathrm{°}`$. Within $`2\mathrm{°}`$ of IRC+10216, little interstellar extinction is detected, since an estimate $`E(BV)0.03`$ mag can be obtained from the maps of Burnstein & Heiles (burstein82 (1982)) and those of Schlegel et al. (schlegel98 (1998)). For such a low extinction, the observations of Guhathakurta & Tyson (guhathakurta89 (1989)) suggest that the diffuse galactic background should be of the order of 27-28 B-mag arcsec-2. This is more than 25 times less than the ISRF surface brightness when it is averaged over the sky (this spatial average is equal to 23.44 $`B`$-mag arcsec-2). Consequently, in order to avoid in our model this fictitious, high background of 23.44 $`B`$-mag arcsec-2, direct photons, which are emitted by the ISRF sphere but are eventually not scattered by the envelope, are ignored.

In the $`UBV`$ range considered here, the model envelope brightness is entirely due to scattered ISRF photons since (i) its optical depth is so large that no photons from the star can escape and (ii) the radiation emitted by the heated dust is insignificant in visible light. Therefore, no central star was implemented in our simulation. For a given set of parameters and for each wavelength, the code returns an image of the externally illuminated envelope. Radial brightness profiles are then derived. To get a profile with a reasonably low noise, around 20 million events (“scattered photons”) are considered per Monte-Carlo simulation and the computation time in parallel mode is about 5 hours on a Compaq AlphaServer equipped with four 525 MHz microprocessors.

## 5 Model results, interpretation and discussion

We consider here the results of modelling when $`standard`$ envelope parameters are adopted (i.e. Table 4, with 0.16 $`\mu `$m grains), and compare them with the observations. We also discuss a model where the grain size is changed to 0.05 $`\mu `$m while keeping the same dust mass loss rate and other parameters. This 0.05-$`\mu `$m model reproduces the observed SED reasonably well (see Fig. 3 of G97), and its radial opacity $`\tau _B`$ is equal to 41.3, which is obtained by scaling the standard opacity with appropriate $`Q_{ext}`$ and $`a`$ values ($`N(r)a^3`$ and $`\tau _BQ_{ext}/a`$, see Eq. 2).

### 5.1 Analysis of the central part of the nebula

The first issue concerns the intensity and shape of the nebula in its central part, within $``$ 50′′ of the center. The $`B`$ band is particularly suitable for consideration because there is no difference between Mathis’s and Mattila’s ISRFs in $`B`$, and because the plateau in $`B`$ is free from a central bright source like in $`V`$. Figure 3 shows that the agreement is poor between the standard model and observations. In order to reproduce the $`B`$-band plateau maximum surface brightness, one has to divide the model intensity by a factor of 2.3 : this factor is clearly too large to be attributed to photometric measurement errors on the plateau brightness. Secondly, concerning the shape itself, the model profile is found to be much broader that the observed one. While the model predicts a radius at half intensity, $`R_{1/2}`$ of 73′′, the observation gives $`R_{1/2}=29`$′′$`\pm `$ 2′′. One notes in Fig. 3 that there is some disagreement between the $`B`$ profiles obtained at OHP and CFH, due to experimental difficulties described in Section 2, but this concerns large offsets and essentially does not affect the $`R_{1/2}`$ value of 29′′. Figure 3 also shows the $`B`$ profile for 0.05-$`\mu `$m grains: in this case, there is a nice agreement on the plateau level, but the model profile is still too broad, with $`R_{1/2}=54`$′′.

The second issue is the wavelength dependence of the nebula shape (shape chromatism). Figure 4 shows the model profiles in $`UBV`$ for the two grain sizes. In each case, we have normalized the $`UBV`$ profiles at the same (arbitrary) maximum value, so that only the shapes are considered. In the standard case (0.16 $`\mu `$m grains), the shape is the same for $`U`$, $`B`$, or $`V`$ and the profiles are exactly superposed. In contrast, if 0.05 $`\mu `$m grains are adopted, there is a strong chromatism in the sense that the $`U`$ profile is much broader than the $`V`$ one. More precisely, for a 0.16 $`\mu `$m grain, $`R_{1/2}=73`$′′ for all filters, whereas for 0.05 $`\mu `$m, one has $`R_{1/2}=71`$′′, 54′′, 29′′ for $`U,B,V`$ respectively. This situation is entirely due to the fact that small grains scatter nearly in the Rayleigh regime, with strong color effects, while larger grains are grey. Because the observations show no chromatism (see Fig. 2) and negligible variation of $`R_{1/2}`$ with wavelength, the larger grain size (0.16 $`\mu `$m) is clearly favoured. It is however true that for small grains the $`V`$-band $`R_{1/2}`$ value fits the observation (29′′) rather well, but there is a large and increasing discrepancy for $`B`$ and $`U`$ bands.

The third issue concerns the colors of the nebula. Figure 5 shows that the observed plateau colors are relatively close to those predicted for the Mattila’s field scattered by 0.16 $`\mu `$m grains. Actually, the agreement is better between the observed colors and the colors of Mattila’s field itself, suggesting pure grey scattering by circumstellar grains, i.e. large sizes. One can also see that if 0.05 $`\mu `$m grains are adopted, none of the fields match the observed colors. Adopting the field of Mathis and 0.16 $`\mu `$m grains is worse than adopting the field of Mattila and the same grains, especially when considering the $`UB`$ index for which the two fields differ. We also note that when one compares the colors of the fields with those of the main sequence stars, the Mathis field seems surprisingly blue in $`UB`$, as blue as an unreddened B8 star. If this was true, the rare stars bluer than B8 would have to compensate the many dwarf or giant stars redder than B8. It is much easier to understand the colors of Mattila’s field, which is located at an average position in $`UB`$ and $`BV`$.

Given the above considerations, it does not appear possible with our data to settle the question of the dominant grain size. Both sizes (0.05 and 0.16 $`\mu `$m) present advantages and disavantages. The small size fits reasonably well the plateau intensity and shape in the V band, but fails for $`B`$ and $`U`$ (shape chromatism) and for the plateau colors. Regarding the large grain size (0.16 $`\mu m`$), the model intensity in $`B`$ is too strong by a factor of 2.3 on the plateau, and also too strong by a factor as large as $`10`$ for $`\theta =70`$′′ (Fig. 3). It is interesting to note that one cannot improve this situation by simply reducing the total opacity (or equivalently the dust mass loss rate) because this produces an even stronger peak at center and disappearance of the central depression, as shown in Fig. 6. We have also explored grains larger than 0.16 $`\mu `$m: with AC1 material, a satisfactory value of $`R_{1/2}`$ is obtained for $`a=0.28`$ $`\mu `$m (then, in the $`B`$ band, $`\tau _{B,\mathrm{rad}}=34.4`$) but the profile does not show any central depression, and the plateau level is still two times too high compared to observation.

A better fit to the profile is obtained by reducing progressively the dust mass loss rate of the standard model in the outer regions. For example, a density decrease in $`r^4`$ at $`r>20`$′′ produces a satisfactory $`R_{1/2}=29`$′′ (a $`r^3`$ law is not enough), although the plateau remains too high.

### 5.2 Analysis of outer regions ($`\theta 70^{\prime \prime }`$)

There is another way of analysing the data: instead of considering the main part of the nebula and the plateau characteristics, one can focus on the outer layers, say at $`\theta =70`$′′. At this radial distance, the surface brightness is certainly less accurately measured than at its center, but its faintness suggests that the envelope can be considered as optically thin at this distance, which greatly simplifies the modelling. Optical thinness is also supported by the fact that the density of background galaxies in a recent very deep $`V`$-band VLT image (de Laverny delaverny02 (2002)) is found to be very similar at $`\theta `$ $``$ 70′′ and farther out. Then, a grain located at this radius is illuminated by almost a pure ISRF for several reasons: i.) the grain sees the central region of the nebula over $``$ 50 degrees FWHM, i.e. over a small solid angle of $``$ 0.6 sr, which is much less than 4$`\pi `$ sr, and ii.) the average surface brightness of the nebula as seen by the grain is fainter than 25.7 mag/arcsec2, which is small against the average ISRF brightness of 23.44 mag/arcsec2 in the $`B`$-band; iii.) most of the directions seen by the grain are optically thin. Consequently, in the optically thin regime and with an isotropic illumination, the surface brightness of the dust envelope at a large angle $`\theta `$ can be written (see, e.g., Guhathakurta & Tyson 1989):

$$S_\lambda (\theta )=\frac{1}{4\pi }f_\lambda \omega \tau _\lambda (\theta ),$$
(3)
where $`S_\lambda (\theta )`$ has the units of ISRF flux $`f_\lambda `$ per steradian, $`\omega `$ is the grain albedo, and $`\tau _\lambda (\theta )`$ is the line of sight optical opacity at offset $`\theta `$. For an unlimited $`r^2`$ envelope, this opacity can be expressed as:

$$\tau _\lambda (\theta )=\frac{3}{16}\frac{Q_e}{a\rho _\mathrm{d}}\frac{\dot{M}_\mathrm{d}}{v_\mathrm{d}d}\frac{1}{\theta },$$
(4)
and finally:

$$S_\lambda (\theta )=\frac{3}{64\pi }f_\lambda \frac{Q_{\mathrm{sca}}}{a}\frac{\dot{M}_\mathrm{d}}{v_\mathrm{d}\rho _\mathrm{d}d}\frac{1}{\theta }.$$
(5)
In principle, the $`\theta ^1`$ dependance of Eq. (5) might be checked on the observed profiles, but in practice, this is hardly possible because of the faintness of the signal and its uncertainties at these large offsets (see Sect. 2, and Figs. 2 & 3).

In Eq. (5), we can reasonably consider that the grain velocity $`v_\mathrm{d}`$, the grain density $`\rho _\mathrm{d}`$, the distance $`d`$ and the ISRF strength $`f_\lambda `$ are relatively well known (Table 3 and 4). Then, because the $`B`$-band observations provide at $`\theta `$=70′′ a surface brightness $`S_B(\theta )`$ of $``$ 4.4 $`10^{20}`$ erg cm-2 s-1 Å-1 arcsec-2 (with an uncertainty of about a factor 2), one gets the following constraint:

$$\frac{\dot{M}_\mathrm{d}}{\mathrm{1.1\; 10}^7\mathrm{M}_{}\mathrm{yr}^1}\times \frac{Q_{sca}}{a_\mu }0.60$$
(6)
For $`a=0.16`$ $`\mu m`$ and $`\dot{M}_\mathrm{d}=\mathrm{1.1\; 10}^7`$ M yr$`^1`$, which is the “standard case”, $`Q_{sca}`$ = 2.16 and $`Q_{sca}/a_\mu =13.5`$. This means that the model flux is found to be $``$ 22 times larger than observed. However, in this case, the line of sight opacity $`\tau _B(70^{\prime \prime })`$ is 0.58, so that the optically thin approximation is invalid here. Fig. 3 illustrates this case more rigourously and indicates a factor of $``$ 9 for the excess of model flux at 70′′.

In order to improve the agreement between model and observations, a lower $`\dot{M}_\mathrm{d}`$ or a lower value for $`Q_{sca}/a_\mu `$ are needed, or both, i.e., less dust or dust with a poorer scattering efficiency. If $`a=0.05\mu m`$, $`Q_{sca}/a_\mu `$ = 4.1 and Eq.(6) implies that $`\dot{M}_\mathrm{d}`$ should be as low as 0.16 $`10^7`$ M yr$`^1`$, and then $`\tau _B(70^{\prime \prime })`$ = 0.042. This “small grain – small dust loss” solution can reproduce the faintness of the nebula at large offsets, but in no case can it account for the central region when extrapolated down to inner layers, because the radial opacity of the whole envelope would be $`\tau _{B,\mathrm{rad}}6`$ ($`\tau _{B,\mathrm{rad}}`$ is equal to $`\frac{1}{\pi }\tau _B(\theta )\theta /\theta _{inn}`$), largely insufficient to produce a plateau. This solution cannot account also for the observed polarized intensity in the near infrared reported by Tamura et al. (tamura88 (1988)).

Larger grains can also be envisaged, and for example with $`a=0.28\mu m`$, $`Q_{sca}/a_\mu `$ is also equal to 4.1; then from Eq.(6), one finds again $`\dot{M}_\mathrm{d}`$ $``$ 0.16 $`10^7`$. In this case, $`\tau _B(70^{\prime \prime })`$ = 0.035 and the envelope radial opacity is found to be $`\tau _{B,\mathrm{rad}}5`$, again largely insufficient to make a central plateau.

### 5.3 Discussion

The above analysis illustrates the fact that it is not easy to reproduce with simple models all the envelope observations. The standard model, which is based on the results of G97 and fits a lot of observations of dust in IRC+10216, especially in the infrared, fails to account for several features of the nebula seen in scattered galactic light. In particular, it seems difficult to simultaneously reproduce the plateau level, its profile, and the faintness of the outer regions.

These difficulties may be due, at least in part, to our simplified treatment of dust. Grain scattering depends primarily on dust optical properties (indices) and the grain size. Scattering also depends on factors like particle porosity and asphericity, and consequently the applicability of the Mie theory. How large the effects of these last two factors are is difficult to estimate without extensive calculations that are beyond the scope of this work. We have also refrained from trying to search for a better agreement between model and observations by simply tuning the optical properties, i.e., modifying at each wavelength the optical indices. The reasons are that G97 used the Rouleau & Martin AC1-type amorphous carbon indices, and, for consistency, the same optical properties should be used here as far as possible; secondly, the indices must obey the Kramers-Kronig relations and cannot be arbitrarily modified in a limited wavelength domain. In summary, our approach has been to consider only two different single particle sizes (0.16 and 0.05 $`\mu m`$), and these two cases correspond to very different albedos and anisotropy $`g`$ factors: compared to the standard case (0.16 $`\mu m`$ grains), the $`0.05\mu m`$ grains have a lower albedo, which gives a very satisfactory fitting of the plateau intensity, but these small grains produce too much chromatism because their albedo strongly depends on wavelength (see above and Appendix C for more details on adopted indices, albedos, etc). It would certainly be interesting to envisage in future work a grain size distribution, although, to our knowledge, no robust observational constraint on circumstellar grain size distributions is presently available (see also the arguments given by G97 for considering only a single grain size).

Our analysis suggests that reducing the dust loss rate at large offsets ($`>25`$ arcsec, corresponding to 1000 yrs ago) is a possible solution to explain the nebula shape and its radial width. Variability of mass-loss with various time scales is a well-known characteristic of AGB stars (see for example the references in Olofsson olofsson99 (1999) or in Marengo et al. marengo01 (2001)), and mass-loss variations for evolved carbon stars is also theoretically predicted (e.g. Wachter et al. wachter02 (2002)). Concerning IRC+10216 itself, observations of its extended envelope in the CO millimeter emission lines have been studied by Huggins et al. (huggins88 (1988)), and more recently by Groenewegen et al. (groenewegen98 (1998)). The last authors conclude, from their analysis of the CO line profiles and maps, that the mass loss is enhanced by a factor of $``$ 5 for $`\theta >50`$′′ compared to inner layers. Such a density enhancement is clearly not seen in our radial profiles, and we checked that if one modifies the standard model by increasing $`\dot{M}_\mathrm{d}`$ by $``$ 5 for $`\theta >50`$′′, a huge bump appears in the envelope profile, and such a bump is definitively not observed (Fig. 2 & 3). Thus, our data tend to suggest that the CO “excess” emission at large offsets is either due, perhaps, to incompletely corrected instrumental lobe effects, or, if real, that they are the consequences of a physical phenomenon different from a density enhancement, as originally proposed by Huggins et al. (huggins88 (1988)).

Concerning the faintness of the nebula, especially at large offsets, it is also interesting to compare the present case of IRC+10216 with the optical studies of high latitude clouds (Guhatakurta & Tyson guhathakurta89 (1989)). These clouds, that absorb and reflect the ISRF, are found to be, in the visible range, about 10 times fainter than predicted by models which reproduce their infrared radiation (with an isotropic ISRF). The optical faintness of these clouds and that of the outer layers of IRC+10216 might be related. Perhaps, if the ISRF is not completely isotropic, both phenomena are due to grain forward scattering of the photons coming preferentially from the galactic plane and escaping away from the observer’s view.

As for the interstellar radiation field, it seems presently premature to draw conclusions about its absolute strength (e.g. in the $`B`$ band) with some accuracy, i.e. an accuracy better than a factor of 2, as was hoped originally when this work started. MR87 proposed that, as a possibility, the field might be a factor $``$ 2.5 smaller than standard, but our calculations with essentially the same envelope characteristics and the same grain radius (0.05 $`\mu m`$) do not confirm this point (Fig. 3). We found no satisfactory reason to explain this inconsistency of results (see Appendix B for more details), but we do find that a fainter field improves the fit with 0.16 $`\mu m`$ or larger grains of AC1-type. Concerning the field colors, our work suggests that the $`UB`$ value of +0.11 found by Mattila (1980a ) is probably more correct than the one (-0.32) from Mathis et al. (mathis83 (1983)).

Clearly, improving the simple model presented here appears desirable in order to reach (hopefully) firmer conclusions. One needs especially to estimate the effects of a grain size distribution, and if possible of grain porosity and asphericity. The effect of a field anisotropy should be investigated. In addition, the influence of envelope homogeneities and structure in discrete shells will have also to be considered, since it is known that clumpiness can significantly modify radiative transfer in envelopes (e.g. Mathis et al. mathis02 (2002)). In parallel with more sophistication in modelling, it would be fruitful to measure the colors of the nebula with less uncertainty, and extend the observations to red wavelengths (when the central variable star is at a minimum) and space ultraviolet. Imaging other similar nebulae might also bring valuable information on the properties of amorphous carbon dust and on the ISRF in the Galaxy.

## 6 Conclusions

We have presented new observations of the IRC+10216 nebula seen in the ambient galactic light, including the first $`U`$ imaging. Additional images in $`B`$ and $`V`$ were obtained as well. These data were used to establish a $`UBV`$ photometric sequence of nearby field stars, and to derive surface brightness profiles of the nebula, and $`UBV`$ colors of the plateau. In blue light, where it is best observed, the nebula has a central depression of 17 %, peaks at $`\theta `$ = 10′′, and its radius at half maximum $`R_{1/2}`$ is 29′′. There is no sign of an increased mass loss at $``$ 50′′, as CO maps may indicate.

The fact that the nebula shape, e.g. its half maximum intensity radius $`R_{1/2}`$, does not depend on wavelength suggests that grains scatter the galactic light in the grey regime. This conclusion is supported by the plateau colors which are close to those of the ISRF as given by Mattila (1980a ). The grain size of 0.16 $`\mu `$m recommended by Groenewegen (groenewegen97 (1997)) is consistent with these aspects of the data; a smaller grain size like 0.05 $`\mu `$m would produce too strong a dependence of $`R_{1/2}`$ on wavelength and too blue a color. However, there remain important discrepancies between the 0.16$`\mu `$m grain model and the observations, essentially concerning the intensity of the plateau and the $`R_{1/2}`$ value, which are better reproduced with 0.05 $`\mu `$m grains.

At first sight, our analysis suggests that the mass loss may be lower in the outer nebula, but reaching firm conclusions is difficult before assessing the effects of a grain size distribution, the envelope structure and perhaps some ISRF anisotropy. Given these uncertainties on grain properties and transfer modelling, there is also no firm evidence that the ISRF incident on IRC+10216 in the visible range is different from that usually taken in the solar neighbourhood.

## Appendix A $`UBV`$ photometry of field stars near IRC+10216

In Table 5, $`UBV`$ photometry is given for 17 objects close to IRC+10216, derived from OHP observations. Coordinates (J2000) are from the APM catalog (Irwin and colleagues; www.ast.cam.ac.uk/apmcat). Objects with flag 1, in the last column, were used, when unsaturated, to calibrate the surface brightnesses of the IRC+10216 nebula in images from OHP, CFH and VLT. Objects with flag 2 are much redder than those with flag 1 ($`BV>1.2`$), and also redder than the nebula; consequently, they were not used for calibration. They are included here for completeness, as is Object 12 with flag 3 which is an elliptical galaxy.

Concerning the photometry, we do not pretend that its quality is the usual one for standard photometric sequences (i.e. $``$ 0.01 mag), and our initial goal was more modestly to reach a $``$ 0.1 mag. uncertainty on colors. The $`UBV`$ data result from aperture photometry made on only 2 to 5 images per color obtained during 3 nights (see the log in Table 1). Calibration on standard fields was limited in accuracy by the weather, which was not quite photometric, except during the night of deep $`U`$ 1-hour exposures. The magnitudes and colors listed in Table 5 were obtained by adopting plausible absorption coefficients for OHP, by including a color term in the transformation from instrumental to standard bandpasses, and by averaging the results from different exposures. The listed uncertainty is the internal error for colors, while it is 0.01-0.02 mag for $`V`$.

A first independent verification of the OHP photometry in $`V`$ and $`U`$ was achieved with the ESO VLT observations. The $`V`$ and $`UV`$ values from VLT were found to be within 0.04 and 0.03 mag of the OHP ones. A second verification was done by considering the color-color diagram of our photometry (Fig. 7). This diagram shows a group of 10 stars with $`BV`$ between 0.45 and 0.8 and located slightly ($``$0.15 mag) above the location of the Population I main sequence (continuous line). This situation is fully consistent with expectations for high latitude (IRC+10216 is at $`b=+44\mathrm{°}`$), faint ($`V`$=14 to 19) stars which statistically probe the metal-weak thick disk.

For instance, Ojha et al. (ojha99 (1999)) provide $`UB`$ vs $`BV`$ diagrams and $`BV`$ histograms for two fields at $`b=+47\mathrm{°}`$, and another exemple of $`UB`$ vs $`BV`$ diagram for $`b=+59\mathrm{°}`$ is given by Yamagata & Yoshii (yamagata92 (1992)). The $`BV`$ edge of the stellar locus from these works is at about 0.4 – 0.45, exactly as in our data, and the $`UB`$ average location is also in fair agreement with our work. Therefore, the final uncertainties of this approximate “sequence” should be probably of $`\pm `$ 0.05 mag. for $`V`$, and $`\pm `$0.05 to 0.1 mag for $`BV`$ and $`UB`$ respectively, at least for (not red) objects flagged 1 used for calibrating the IRC+10216 surface brightness. The sequence is not the main source of error in the nebula surface brightness photometry, whose accuracy (0.1 mag in $`V`$, 0.2-0.25 mag in color indices) is more limited by its intrinsic faintness and sky level subtraction.

## Appendix B Quantitative verifications of the numerical code

Some numerical tests have been done in order to verify the consistency of the results returned by the code. They are in complete agreement with other computations or analytical relations.

First, the computed SED of IRC+10216 with the parameters of Table 4 is in excellent agreement with that published by G97. This increases confidence in the resolution of the radiative equilibrium equations.

Secondly, in order to check the computation of azimuthally-averaged brightness profiles (i.e. simulations of the ISRF illumination and scattering in the envelope), we considered the ideal case of an optically thin homogeneous sphere illuminated by an isotropic ISRF. We assume that the dust density is constant in the sphere of external angular radius $`\theta _{\mathrm{out}}`$. Because of optical thinness, each grain is illuminated by the ISRF itself, and it can easily be shown that the surface brightness of such a sphere at an angle $`\theta `$ from its center can be expressed as:

$$S_\lambda (\theta )=\frac{1}{2\pi }f(\lambda )\frac{Q_{\mathrm{sca}}(\lambda )}{Q_{\mathrm{ext}}(\lambda )}\tau _{\lambda ,\mathrm{rad}}\sqrt{1(\frac{\theta }{\theta _{\mathrm{out}}})^2}$$
(7)
where $`Q_{\mathrm{sca}}(\lambda )`$ and $`Q_{\mathrm{ext}}(\lambda )`$ are respectively the scattering and extinction efficiencies at wavelength $`\lambda `$, $`f(\lambda )`$ the ISRF flux and $`\tau _{\lambda ,\mathrm{rad}}`$ the radial optical thickness of the sphere (i.e. from center to outer radius). The brightness $`S_\lambda (\theta )`$ has the units of $`f(\lambda )`$ per steradian. An agreement better than 1% is obtained between the numerical and analytic brightness profiles for $`\tau _{\lambda ,\mathrm{rad}}0.01`$, i.e. when the single scattering assumption remains correct. This test validates our treatement of ISRF illumination with an external sphere.

Thirdly, for an optically thin circumstellar envelope with a radial density varying as $`r^2`$ , the brightness profile can also be analytically estimated (for $`\theta \theta _{\mathrm{inn}}`$):

$$S_\lambda (\theta )=\frac{1}{4}f(\lambda )\frac{Q_{\mathrm{sca}}(\lambda )}{Q_{\mathrm{ext}}(\lambda )}\tau _{\lambda ,\mathrm{rad}}\frac{\theta _{\mathrm{inn}}}{\theta }q\text{}$$
(8)
where $`q^{}=1`$ if no outer limit is assumed:

$$q^{}=(1\frac{R_{\mathrm{inn}}}{R_{\mathrm{out}}})^1\frac{2}{\pi }arctan(\frac{\sqrt{1u^2}}{u})$$
(9)
with $`R_{\mathrm{inn}}`$ and $`R_{\mathrm{out}}`$ being the inner and outer envelope radius respectively, $`\theta _{\mathrm{out}}=R_{\mathrm{out}}/d`$ and $`u=\theta /\theta _{\mathrm{out}}`$. In this case also, the code returns very good results in comparison to this analytical brightness profile, provided that $`\tau _{\lambda ,\mathrm{rad}}`$ is low enough.

We have also examined the case of similar but optically thicker spheres. For homogeneous spheres, our simulations show that the brightness profiles depart from the above equation (C.1) when multiple scattering plays a role, i.e. around $`\tau _{\lambda ,\mathrm{rad}}>0.1`$. An almost constant brightness profile with $`\theta `$ is obtained for $`\tau _{\lambda ,\mathrm{rad}}1`$ and bright-rimmed spheres with dark cores are found for optically thicker cases. We also found that for optically thick spheres with radial density varying as $`r^2`$, the surface brightness shows a minimum at the center, reaches a maximum at some angular distance from the center and then decreases strongly. All these results are entirely consistent with previous similar simulations such as those found in Witt & Stephens (witt74 (1974), see their Fig. 1 & 2).

Finally, we have attempted to exactly reproduce the numerical results of Martin & Rogers (martin87 (1987)). One can first note that their Fig. 8 shows a significant chromatism of the plateau profile between the $`I`$ and the $`V`$ bands (0.81 $`\mu `$ and 0.55 $`\mu `$) that is qualitatively consistent with our findings explained in Sect. 5. At 0.55 $`\mu m`$, for a standard ISRF, their modelled plateau is at 2.7 fu (flux unit, where 1 fu is $`10^{22}`$ W m-2 Hz-1 sr-1), extends to $`\theta 5`$′′, and has $`R_{1/2}`$ 25′′.

Their model has essentially the same physical and geometrical properties as ours, but unfortunately their optical grain properties are not given in an explicit way. It is mentioned in their Sect. III-i that in $`V`$, $`\omega `$=0.38 and $`g`$=0.10. In addition, from their Fig. 1 where $`Q_{abs}/a`$ and $`Q_{sca}/a`$ are plotted, one can estimate $`Q_{abs}0.21`$ and $`Q_{sca}0.15`$. Using the Mie theory, with $`a`$=0.05$`\mu `$m and $`\lambda `$=0.55$`\mu `$m, one finds corresponding indices $`n`$=2.470 $`k`$=0.264, for which $`\omega `$=0.384 and $`g`$=0.15, $`Q_a`$=0.221, $`Q_s`$=0.138, $`Q_e`$=0.359, so that only $`g`$ is discrepant (we have found no reason for that discrepancy). The model of Martin & Rogers has also $`\tau _{11\mu }`$=0.8 and from their fig. 1, $`Q_e`$(11$`\mu `$)=$`Q_a`$=0.0125. Therefore, $`\tau _{0.55\mu }`$ = $`\tau _{11\mu }`$ $`\times `$ $`Q_e`$($`V`$)/$`Q_e`$(11$`\mu `$) = 23.0, in fair agreement with the value of 21 mentioned by them.

With $`\tau _{0.55\mu }`$=23 and the $`n`$,$`k`$ values found above, our code finds for a standard ISRF a very peaked central maximum at 2.75 fu, with no plateau and $`R_{1/2}`$ $``$ 15′′. Only when $`\tau _{0.55\mu }`$ is increased to 30 does one get a plateau with maximum intensity of 2.0 fu, extended over 5′′ and with $`R_{1/2}`$ 36′′. These two cases bracket the result of MR87, but are significantly different. Therefore, there are some differences between the results of Martin & Rogers and those obtained with our code, but it is difficult to say whether they are due to differences in the codes themselves or due to differences in indices.

## Appendix C Grain optical properties

Table 6 lists the optical properties adopted in this work (amorphous carbon of type AC1 from Rouleau & Martin rouleau91 (1991)).

Acknowledgments: The authors thank the anonymous referee and P.J. Huggins for comments that helped us improve the paper. We also thank B. Gladman for careful reading of the manuscript. Part of this work has been performed using the computing facilities provided by the program “Simulations Interactives et Visualisation en Astronomie et Mécanique (SIVAM)” at Observatoire de la Côte d’Azur. We acknowledge support from the CNRS Program Physico-chimie du milieu interstellaire (to N.M.), the CNRS-INSU Actions thématiques innovantes and the Ministère de l’Education Nationale et de la Recherche (to P.d.L.).