text:
# The mineral composition and spatial distribution of the dust ejecta of NGC 6302 Based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA Based on observations made with ESO Telescopes at the La Silla or Paranal Observatories under programme ID 67.D-0132(A)

## 1 Introduction

All low and intermediate mass stars end their life on the Asymptotic Giant Branch (AGB) by ejecting their entire H-rich envelope through a slowly expanding, dusty wind. After the AGB, the central star quickly increases its effective temperature to values high enough to begin ionizing its AGB ejecta: a planetary nebula (PN) is born. Depending on the mass and thus luminosity of the star, the transition from AGB to PN can go very fast: for the most massive objects, in less than 1000 years ionization of the nebula begins. These massive objects therefore are characterized by dense AGB remnants and a very hot luminous central star.

NGC~6302 is probably one of the best studied PNe with a massive progenitor. A recent determination of the mass of the ionized nebula is about 2 $`M_{}`$ (Pottasch & Beintema 1999), based on a distance determination of 1.6 kpc (Gómez et al. 1993). Dust mass estimates also indicate a very massive shell ($`M_\mathrm{d}=0.02M_{}`$), using a distance of 2.2 kpc (Gómez et al. 1989). These distance determinations rely on VLA observations of the expansion of the nebula over a 2.75 yr period, assuming an outflow velocity of 13 km s-1. In principle, this is a very reliable method, however the epoch over which the nebula is observed, is very short. Therefore the increase in size is very small and hard to measure, indicating that the error bars on these distance determinations are still very large. Instead, we adopt a distance of 0.91 kpc, based on emission-line photometry of NGC~6302, from which the $`B`$ and $`V`$ magnitudes, the luminosity and the distance can be derived (Shaw & Kaler 1989). The large number of luminous PNe studied gives confidence in the distance determination from this statistical method.

The nebular abundances indicate that NGC~6302 is a type I nebula, consistent with a massive progenitor. Casassus et al. (2000) estimate that the progenitor mass is 4–5 $`M_{}`$. The morphology of the nebula observed at optical wavelengths is highly bipolar, pointing to non-spherical mass loss on the AGB, resulting in a dusty torus in the equatorial region (Lester & Dinerstein 1984). The inclination angle of the system is $`45^\mathrm{o}`$ with respect to the line-of-sight (Bohigas 1994). The temperature of the central star is very high: Casassus et al. (2000) mention a temperature of $`\mathrm{250\hspace{0.17em}000}`$ K, while Pottasch et al. (1996) arrive a temperature of $`\mathrm{380\hspace{0.17em}000}`$ K. Although there is some uncertainty about the distance and therefore masses and luminosities involved, everything points to a rather massive and luminous progenitor.

The Infrared Space Observatory (ISO) spectrum has been presented in several papers (Barlow 1998; Molster et al. 2001a, 2002b) and is characterized by a wealth of narrow solid state features in the 20–70 $`\mu `$m spectral range caused by circumstellar dust in the AGB remnant, as well as by strong emission lines from a multitude of fine-structure lines originating from the ionized gas in the nebula. The dust bands have been identified with crystalline silicates and a number of other components (Koike et al. 2000; Molster et al. 2001a, 2002a). Kemper et al. (2002b, hereafter referred to as Paper I) have reported on the detection of carbonates in the dust shell, based on the identification of broad features at $``$60 and $``$92 $`\mu `$m.

Optical and near-infrared images have already shown that the dust distribution around NGC~6302 is rather complex (Lester & Dinerstein 1984; Schwarz et al. 1992). The ISO spectroscopy supports this, because a broad dust temperature range is required to explain the shape of the spectral energy distribution. In addition the dust composition seems to be very complex, as there is evidence for a mixed chemistry by the presence of both oxygen-rich dust and carbon-rich dust features, the latter in the form Polycyclic Aromatic Hydrocarbons (PAHs) (see Molster et al. (2001a) for a discussion on the origin of this dichotomy).

In order to reconstruct the mass loss history of NGC~6302, including the geometry and composition of the AGB wind, infrared spectroscopy and imaging are needed. Unfortunately, the ISO data lack spatial information, limiting the analysis to the bulk dust composition. Mid-infrared imaging can reveal the present-day geometry of the dust envelope, which puts limits on the mass loss history.

This paper is organized as follows: In Sect. 2 the ground-based mid-infrared imaging and spectroscopy is presented. The ISO spectroscopy is discussed in Sect. 3, which describes our laboratory data of carbonates (Sect. 3.1) and a model fit to the observed spectrum (Sect. 3.2). In Sect. 4 we propose a possible geometry of NGC~6302 and discuss the astronomical relevance of carbonates. Our results are summarized in Sect. 5.

## 2 The TIMMI2 observations

### 2.1 N- and Q-band imaging

We have observed NGC~6302 using the TIMMI2 mid-infrared imaging spectrograph attached to the 3.6m telescope at the European Southern Observatory, La Silla, Chile. For a description of the instrument see Reimann et al. (2000). The observations were carried out on the night of June 17/18, 2001. Observing conditions were not optimal, with occasional clouds and variable seeing (between about 0.6′′ and 1.4′′). Imaging was done in four photometric bands, centred at 9.8 (N2), 10.4, 11.9, and 20 (Q) $`\mu `$m. In addition, we obtained two 8.35-11.55 $`\mu `$m long-slit spectra. We observed HD~81797 and HD~169916 as point sources in order to determine the shape of the point spread function, as well as for flux calibration purposes, which is only performed in case of spectroscopy.

To correct for background, chopping and nodding has been performed. Apart from the Q-band image, the chopping and nodding images fell all inside the frame. For the various N-band images we used a lens scale corresponding to a pixel size of 0.3′′/pixel, whereas the Q-band image was obtained with a pixel size of 0.2′′/pixel. The dimensions of the detector are 320$`\times `$240 pixels of which the central 300$`\times `$220 pixels were used for further analysis. The chop throw used for the images in the N-band is 30′′, with a nodding offset of 30′′ perpendicular on the chop direction. For the Q-band image, we applied a chop throw of 40′′, and a nodding offset of 40′′ in the same direction as the chop.

Both the positive and negative images were used for the resulting image. The final image has been sharpened by means of a deconvolution using the point spread function (PSF) derived from HD~81797 for the Q-band. Unfortunately, for the N2, 10.4 and 11.9 filters no standard star was measured in the same night with the same settings, and we used observations of standards of another night. For the N2 and 10.4 filters no standards were available at all. Therefore, we used HD~169916 measured with the 11.9 $`\mu `$m filter for all three filters for the PSF, assuming that the PSF would be similar in all three N-bands. A comparison with other N2-band images (with different settings, but the positive and negative images did overlap each other) showed that our assumption is not unreasonable.

Fig. 1 shows the final TIMMI2 images at N2, N10.4, N11.9, and Q band, as well as the ”raw” (i.e. before deconvolution) N2 image, and an optical VLT image of the nebula for comparison (ESO 1998). The ”raw” N2 image shows a bright elongated region of approximately 10′′ by 4′′, and a much fainter region offset to the south and west. There is also evidence for faint arcs of emission that are reminiscent of the X-shaped nebula seen in the optical image. In fact, the mid-IR emission has a morphology which is strikingly similar to that seen in the optical (see Fig. 2). This is surprising since the mid-IR emission originates from (warm) dust, while the optical emission is dominated by ionised gas. The equatorial torus, seen as a waist in the optical images, must contain very cold dust which does not yet contribute to the flux at 10 $`\mu `$m. We will return to this point when we discuss the ISO spectrum.

The deconvolved N2, N10.4 and 11.9 $`\mu `$m images in Fig. 1 are very similar and show a rather ”blobby” appearance. Note that the faint structure in the background is a common artifact of image deconvolution. However, the structure in the central area is real, since we see similar substructure in different observations deconvolved with different point spread functions. Also the 6cm radio map of NGC 6302 shows substructure in this region (Gómez et al. 1989). Our images are very similar to the 3.3 $`\mu `$m image taken by Casassus et al. (2000), taking into account that the orientation of their image is accidentally rotated by 180o (Casassus, priv. comm.).

The deconvolved Q band image (Fig. 1) tends to show more spherical emission, although the eastern part of the nebula is still dominating the emission. Due to the reduced spatial resolution at 20 $`\mu `$m it is not possible to say if the emission is equally blobby as seen at 10 $`\mu `$m.

### 2.2 N-band spectroscopy

In Fig. 3 we show the spectra taken at the slit positions indicated in Fig. 4. We used HD~169916 as a reference spectrum, assuming that this star has a $`F(\nu )\nu ^2`$ flux distribution. The comparison with other calibration observations showed that the absolute flux calibration is only accurate up to about 40%. The relative flux calibration is much more accurate. The wavelength calibration provided by ESO on the TIMMI2 home page appeared to be inaccurate, and we determined a new wavelength calibration, based on the \[Ariii\] and \[Siv\] line, assuming a linear pixel to wavelength dependence. The fact that the \[Navi\] line and the PAH features are now at the right wavelengths shows that our assumption is valid (using the calibration files provided by ESO the 11.3 $`\mu `$m PAH feature was found at 11.7 $`\mu `$m). The comparison with the ISO spectrum (dashed line in Fig. 3), gives confidence in the new wavelength calibration.

The spectrum is dominated by the forbidden emission lines of \[Siv\] and \[Ariii\]. However the PAH features at 11.3 and also 8.6 $`\mu `$m are also found. The noise around 9.5 $`\mu `$m is due to telluric lines.

The east and west slit are both probing the eastern lobe of the nebula (Fig. 4). The difference spectrum is very flat, apart from the \[Siv\] line. The fact that even around 8.6 and 11.3 $`\mu `$m it is flat, indicates that the distribution of the PAHs is quite constant in the bright eastern lobe. This is in agreement with the 3.3 $`\mu `$m image obtained by Casassus et al. (2000).

## 3 Analysis of the ISO spectrum

The ISO SWS and LWS spectra have been presented by Barlow (1998) and Molster et al. (2001a, 2002b) and in Paper I. We briefly repeat the overall characteristics of the spectrum. The SWS and LWS apertures, which have sizes in the range 14$`\times `$20′′ to 80′′, cover a large fraction of the nebula and therefore provide no detailed spatial information. The spectrum, shown in Fig. 5, is dominated by a strongly rising spectrum which peaks near 55 $`\mu `$m, and numerous, sometimes rather narrow emission bands. Both the continuum and these emission bands are caused by the circumstellar dust. In addition, strong, spectrally unresolved emission from fine-structure lines are evident, which are caused by the ionized nebular gas. These lines can be used to determine the chemical abundance of the gas and properties of the central star (e.g. Pottasch et al. 1996; Pottasch & Beintema 1999). In this work we focus on the analysis of the dust component.

### 3.1 Laboratory spectroscopy of carbonates

Molster et al. (2001a) provide a useful inventory of the solid state emission bands in the ISO spectrum of NGC~6302, and we will not repeat their analysis here. While many emission bands were successfully allocated to (mostly oxygen-rich) dust components, several prominent emission bands, notably near 29, 48, 60 and 90 $`\mu `$m remained partially unidentified. We have investigated possible carriers for the $``$92 $`\mu `$m band, reported to be present in NGC~6302 and NGC~6537 (Molster et al. 2001a, 2002b). During the course of our laboratory studies of cosmic dust analogues, we noted that the carbonate calcite (CaCO3) has a strong resonance at $``$92 $`\mu `$m. This prompted us to carry out a systematic laboratory study of transmission spectra of several carbonate species. A series of carbonates of possible cosmic relevance has been chosen for laboratory investigations, including the natural minerals ankerite (origin unknown), aragonite (Bohemia), dolomite (Eugui, Navarra, Spain), calcite (origin unknown), magnesite (Steiermark, Austria) and siderite (Schönbrunn, Vogtland, Germany). The purity and homogeneity of the natural samples are investigated by scanning electron microscopy (SEM) and energy dispersive X-ray analysis (EDX). The macroscopic phase homogeneity of the carbonate samples was measured by X-ray diffraction. The chemical composition and the content of impurities are presented in Table 1.

The minerals dolomite, calcite, aragonite and magnesite turned out to be very pure and homogeneous. The impurities in siderite and ankerite are the corresponding carbonates mentioned in Table 1. These are homogeneously distributed within the main mineral component.

Small carbonate grains with an average size smaller than 2 $`\mu `$m have been produced by grinding and sedimentation in water-free acetone. In order to measure the infrared spectrum over a wide spectral range the small grains were embedded in transparent material like KBr and Polyethylene (PE). Infrared spectroscopy was performed using a Bruker 113v Fourier transform spectrometer in the wavelength range between 2-200 $`\mu `$m. The infrared spectral behaviour of the carbonate minerals is demonstrated in Fig. 6, and can be found in digitized form in the Jena-St. Petersburg Database of Optical Constants (Henning et al. 1999).

Magnesite, siderite, ankerite, dolomite, and calcite belong to the carbonates of the calcite group and crystallize in the trigonal crystal system. Carbonates containing cations with larger radii than Ca2+ crystallize in the rhombic crystal system. CaCO3 is dimorph and can exist in both structures as calcite and aragonite. The typical crystalline structure of the carbonates in the calcite group is characterized by an alternating sequence of cation and anion layers. In both magnesite and calcite, the plane carbonate group consists of three oxygen located at the corners of an equilateral triangle and a carbon in the centre. The bonding between C and O is covalent whereas the bonding between the cations and the carbonate anion (CO32-) is mainly ionic. Each cation is octahedrally coordinated by six O atoms of the adjacent CO32- ion. The crystals of the calcite group with mixed metal composition – like dolomite and ankerite – show an alternating arrangement of the calcium and magnesium or calcium and iron layers. In contrast to the calcite structure the Ca2+ ion in aragonite is coordinated by 9 O-atoms.

The infrared optical properties are determined by the crystalline structure of the minerals. The spectral behaviour of carbonates has been measured several times in previous studies (e.g. Hellwege et al. 1970; Penman 1976; Orofino et al. 1997; Yamagishi et al. 2001). The measured mid-infrared bands arise mainly from the carbonate anion. The planar structure of this anion and its C3 symmetry gives rise to three infrared active modes at about 7, 11 and 14 $`\mu `$m. The 7 $`\mu `$m band can be attributed to the C=O stretching mode whereas the 11 and 14 $`\mu `$m bands are caused by the out-of-plane and in-plane bending modes. The low frequency bands beyond 25 $`\mu `$m are due to translation of the metal cations and represent motions perpendicular and parallel to the plane of the carbonate anions (Hellwege et al. 1970).

### 3.2 Dust model fit

#### 3.2.1 Description of the model

We use the identified solid state components to quantitatively fit the entire ISO spectrum. Inspection of the ISO data shows that all dust bands are in emission, suggesting that the dust shell is optically thin at infrared wavelengths. This substantially simplifies the analysis. Given the complex geometry of the nebula, the multiple dust components present in the spectrum, and the still incomplete knowledge of the dust species present, we decided to focus our analysis on determining the relative mass contributions of the known dust species, rather than attempting a full 2-dimensional radiative transfer calculation. We use the optically thin nature of the dust at mid- and far-infrared wavelengths to constrain the dust mass of the components using mass absorption coefficients. The thus derived dust masses will not depend on the assumed geometry as long as the emission is optically thin. The problem then reduces to deriving the temperature over mass distribution of the various dust components (e.g. Bouwman et al. 2000). The advantage of using mass absorption coefficients in an optically thin dust shell, is that it does not require knowledge of the location of the dust.

We derive the dust mass and temperature distribution under the following assumptions: i) the dust shell is optically thin at mid- and far-infrared wavelengths; ii) all grains have the same size, i.e. 0.1 $`\mu `$m; iii) all grains are spherical; iv) all grains are of homogeneous composition; v) the number density of dust grains as a function of distance is a simple power law (see Eq. (4), where $`\rho _0=\frac{4}{3}\pi a^2\rho _\mathrm{d}n_0`$); vi) the dust species all have the same temperature ranges, limited by the temperatures at the inner radius $`T_0`$ and at the outer radius $`T_{\mathrm{max}}`$, which are free parameters. The result of our analysis is a mass over temperature distribution $`M_\mathrm{d}(T)`$ of all identified dust species. In Sect. 4.1.1 we will discuss the validity and the impact of these assumptions.

In the general case of an optically thin dust cloud, which specifically applies to an optically thin dust shell, the specific intensity emitted by the dust in the cloud can be written as

$$I_\nu (T)=_V\alpha _\nu B_\nu (T)𝑑V\text{for}\tau 1$$
(1)
where $`\alpha _\nu `$ is the specific absorptivity. By integrating over volume $`V`$ this can be written as

$$I_\nu (T)=N\sigma _\nu B_\nu (T)\text{or}I_\nu (T)=M_\mathrm{d}(T)\kappa _\nu B_\nu (T)$$
(2)
with the absorption cross-section can be described in terms of the geometric cross-section and absorption efficiency, according to $`\sigma _\nu =\sigma _{\mathrm{geom}}Q_\nu `$. Following for instance Bouwman et al. (2000), $`M_\mathrm{d}(T)`$ is the dust mass at temperature $`T`$. The mass absorption coefficient is given by $`\kappa _\nu `$ and the number of grains in the volume is represented by $`N`$.

When we realize that we can write for the specific flux $`F_\nu =I_\nu /D^2`$, and under the assumption that the dust grains are identical spherical grains with a radius $`a`$, we find

$$F_\nu (T)=\frac{1}{D^2}\frac{1}{\rho _\mathrm{d}}\frac{3}{4}\frac{1}{a}M_\mathrm{d}(T)Q_\nu B_\nu (T)$$
(3)
with $`D`$ the distance to the dust cloud and $`\rho _\mathrm{d}`$ the density of the material in the dust grain.

In case of a circumstellar dust shell, we may assume that the temperature and density distribution follow a power law:

$$T(r)=T_0\left(\frac{r}{r_0}\right)^q\rho (r)=\rho _0\left(\frac{r}{r_0}\right)^p$$
(4)
where $`r`$ represents the distance to the central star, and $`r_0`$ the inner radius of the dust shell. Using $`dM_\mathrm{d}=4\pi r^2\rho (r)dr`$, we can express $`M_\mathrm{d}`$ in terms of the temperature distribution in the dust shell, according to

$$M_\mathrm{d}(T(r))=\frac{4\pi \rho _0r_{0}^{}{}_{}{}^{3}}{3p}\left(\frac{T(r)}{T_0}\right)^{\frac{3p}{q}}$$
(5)
For spherical dust grains the density at the inner radius is $`\rho _0=\frac{4}{3}\pi a^3\rho _\mathrm{d}n_0`$, with $`n_0`$ the dust grain number density at the inner radius, which is different for each identified dust component. The density $`\rho _\mathrm{d}`$ corresponds to the density of the grain material itself. This can be combined with Eqs. (5) and (3) to obtain an expression of the specific flux:

$$F_\nu (T)=\frac{1}{D^2}\frac{4\pi a^2r_{0}^{}{}_{}{}^{3}n_0}{3p}\left(\frac{T(r)}{T_0}\right)^{\frac{3p}{q}}Q_\nu B_\nu (T)$$
(6)
where the observable total specific flux $`_\nu `$ is given by

$$_\nu =_{T_0}^{T_{\mathrm{max}}}F_\nu (T)𝑑T$$
(7)
which can be compared directly to the ISO spectroscopy.

#### 3.2.2 Fit parameters and dust composition

Following the analysis of Molster et al. (2001a), we assumed that enstatite (MgSiO3), forsterite (Mg2SiO4) and crystalline water ice are present, in addition to amorphous silicates. Diopside ((Ca,Mg)SiO3) is believed to contribute to the 60 $`\mu `$m complex (Koike et al. 2000). Diopside is in the form of a solid solution between wollastonite (CaSiO3) and enstatite (MgSiO3) in the number ratio 46:54, i.e. Wo46En54. In addition, we considered metallic iron as a dust component, which is found to be present in the dust shells of evolved stars (Kemper et al. 2002a), and the carbonate species presented in this work. Under the assumption that all dust grains are spherical and have a radius $`a`$ of 0.1 $`\mu `$m, we have derived $`Q_\nu `$ values from laboratory measurements of amorphous olivine (Dorschner et al. 1995); metallic iron (Henning & Stognienko 1996); forsterite and clino-enstatite (Koike et al. 1999); diopside (Koike et al. 2000); water ice (Bertie et al. 1969; Warren 1984) and carbonates (this work).

Using Eqs. (6) and (7) the far-infrared spectrum has been fitted, under the assumption that all dust species are found in the same temperature regimes. Thus, the free parameters to be fitted are the number density at the inner radius $`n_0`$, which differs for each dust component and scales with the mass contained in that component; the temperature range defined by $`T_0`$ and $`T_{\mathrm{max}}`$; and the slopes $`p`$ and $`q`$ of the density and temperature power laws. For the distance $`D`$ we have used the measured value of 0.91 kpc (Shaw & Kaler 1989) and the grain size $`a`$ is assumed to be 0.1 $`\mu `$m. The inner radius $`r_0`$ follows from the values of $`T_0`$ and $`q`$, and is therefore not a fit parameter. If $`p`$ and $`q`$ are both chosen to be $`1/2`$ the best results are obtained. Inspection of the spectrum shows that the temperature distribution of the dust must be very wide, which is incompatible with a simple $`r^2`$ density distribution of the dust. We attempted to fit the spectrum with a single range of temperatures but were not successful; apart from a cool dust component (30–60 K), peaking near 55 $`\mu `$m, we had to introduce a second dust component at 100–118 K dominating the spectrum between 10 and 30 $`\mu `$m (see Fig. 7). The results of the fits are shown in Table 2, where the mass of each dust component for both temperature ranges is given. These masses are constrained by $`n_0`$ which is separately fitted for all identified dust species and both temperature ranges, and are obtained by integrating Eq. (5) over the considered temperature range.

The warm dust component could not consist of amorphous olivine, because the strong resonances due to the Si-O stretching and O-Si-O bending modes at 9.7 an 18 $`\mu `$m are not observed. Iron is used as the carrier of the warm continuum, but as PAHs have been detected in NGC~6302 (Molster et al. 2001a), amorphous carbon, instead of iron, could also explain the continuum emission. A mixed chemistry is not uncommon for evolved stars, it is for example seen in IRAS~09425$-$6040 (Molster et al. 2001b). We conclude that the composition of the warm dust is different from the composition of the cool dust.

The exact composition of the cold amorphous dust is unknown. Amorphous olivine is used to produce the required continuum emission, but other abundant species with a smooth absorption spectrum at far-infrared wavelengths, such as pyroxene and iron, could be contributing as well. The mass fractions of forsterite and enstatite are lower limits. We used the laboratory data leading to the best spectral match with the observations, i.e. those of Koike et al. (1999).

The previously unidentified 92 $`\mu `$m band (Molster et al. 2001a) can be fitted very accurately with calcite (CaCO3, trigonal), see Fig. 8. The strong resonances of calcite at shorter wavelengths (see Fig. 6) are not discernible in the ISO spectrum, indicating that it is not a significant component in the warm dust. Aragonite (CaCO3, rhombic) also has a resonance at $``$92 $`\mu `$m, and may therefore contribute to the observed feature in NGC~6302. However it is not very likely to exist in a Planetary Nebula. Although it has the same chemical composition as calcite, its lattice structure is different, and it can only exist under relatively high pressure and temperature.

The 60 $`\mu `$m feature can be accurately fitted using absorption efficiencies of dolomite, diopside and water ice (see Fig. 9). Ankerite (FeCO3) has a resonance at approximately the same wavelength (see Fig. 6) and could in principle contribute to the 60 $`\mu `$m complex. However, solar system carbonates are almost exclusively in the form of dolomite and calcite (Johnson & Prinz 1993; Brearley & Jones 1998) and therefore, we concentrated our analysis on dolomite.

The results of our model fit are shown in Fig. 8 and Table 2. Since the dust shell is optically thin, the relative mass fractions can be derived from the model fit. The absolute dust mass of each mineralogical component scales with the distance, and if 910 pc is adopted for the distance (Shaw & Kaler 1989), we arrive at the dust masses as indicated in Table 2. We derived that the depletion of calcium into the identified calcium bearing species – calcite, dolomite and diopside – is about 30% (Paper I), if we assume a dust/gas mass ratio of 1/100 and adopt a calcium abundance of $`[`$Ca/H$`]=2.210^6`$ (Snow 1984). From the fit to the 60 $`\mu `$m complex the fraction of water contained in the solid phase can be calculated. Adopting a mass ratio of $`110^3`$ between H2O and H2 (González-Alfonso & Cernicharo 1999) we find that $``$ 10% of the water is contained in water ice. This is consistent with the results of Hoogzaad et al. (2002). The presence of OH molecules (Payne et al. 1988) indicates that water should also be present in the vapour phase.

## 4 Discussion

### 4.1 Geometry and composition of the circumstellar dust shell

In Sect. 3.2.2 we derived that the dust in the circumstellar environment of NGC~6302 consists of two temperature components. The cold component (30–60 K) has a mass of 0.050 $`M_{}`$ and contains silicates, water ice and carbonates, indicative of a oxygen-rich chemistry. The warm component (100–118 K) contains only $`1.310^4`$ $`M_{}`$, which is mainly caused by a featureless species like carbon or iron. About 5% of the warm component is in the form of silicates. The large difference in chemical composition between the warm and the cold component and the lack of continuity in the temperature distribution suggest that the two components have a different spatial distribution. This then allows a mixed chemistry; regions of C-rich dust could exist in a PN with a strong silicate spectroscopic signature. The presence of PAHs in NGC~6302 (first reported as UIR bands by Roche & Aitken 1986) supports the idea of a mixed chemistry, both the warm, diffuse component and the PAHs could reside in the same regions of the nebula.

The large mass of the cold component suggests that it is present in a circumstellar torus, which leaves the polar regions for the warm carbon-rich component. Optical images of NGC~6302 show that the eastern lobe is brighter than the western lobe (see Fig. 1 and also Schwarz et al. 1992), suggesting that the western lobe is partly obscured by the circumstellar torus. This leads to a geometry as indicated in Fig. 10. The dusty torus around the central star is seen under an inclination angle $`45^\mathrm{o}`$, and partially obscures the western lobe. In return, the eastern part of the torus is located behind the eastern outflow. The dusty torus contains a high column density of dust towards the central star, and therefore, most of the dust is effectively shielded from heating by the central star and remains relatively cool. The dust in the polar regions however, is directly radiated by the central star and thus has a higher temperature.

Our observations, as well as previously published work, is in agreement with the proposed geometry. The torus is optically thick in the UV, optical and near-infrared. Therefore, the optical images only show the scattered light from the central star in the polar regions. We have convolved the warm and cool dust spectra modelled in Sect. 3.2.2 with the TIMMI2 N11.9- and Q-band transmission profile and corrected for the atmospheric transmission in order to determine the contribution of the dust components in both bands. We find that in the N11.9-band, only 0.2% of the flux comes from the cold component, and consequently 99.8% originates from the warm dust component. For the Q-band these numbers are 29.0% and 71.0% respectively. Since the peak in intensity in the N-band images and the optical image coincide (see Fig. 2) we conclude that indeed the warm dust is located in the polar regions. The Q-band imaging reveals a different picture; one third of the flux originates from the cool dust component, and two third from the warm dust. Both the warm dust in the eastern lobe, and the cool dust in the western part of the torus contribute to the observed image, leading to a more or less spherical structure in this wavelength band (Fig. 1).

N-band spectroscopy of the eastern lobe shows PAH emission. From Fig. 6 of (Casassus et al. 2000, figure is rotated by 180o) it is clear that the PAH emission intensity peak coincides with the intensity peak in the optical and N-band, indicating that the PAH emission also arises from the lobes, where UV radiation from the central star can penetrate. The UV radiation field also explains the ionization observed in the lobes (Fig. 3). Persi et al. (1999) also notice that the emission of warm dust grains and forbidden line transitions emerges from the same region. The 6cm data presented by Gómez et al. (1989) traces the location Hii region. In their Fig. 3 the region confined by the torus is clearly visible.

On the other hand, the optically thick torus around the central star effectively shields the molecular material in the torus from the high energy radiation. The infrared radiation field in the torus will pump OH molecules present and cause the formation of OH masers. Indeed, OH maser emission is observed in NGC~6302, on the west side of the object (Payne et al. 1988), indicating that it arises from part of the torus located at the front side of the object. This is consistent with the fact that predominantly the blue-shifted part of the line profile is seen (Payne et al. 1988).

#### 4.1.1 Validity of model assumptions

The geometry described here is consistent with the assumptions we have used to fit the ISO spectrum, and with the TIMMI2 images. Here we will discuss the validity of the assumptions used in this model.

i) We assumed that the dust is optically thin at large wavelengths. The TIMMI2 images indeed suggest that this is true: the torus is still visible as a dark lane in the N-band images but starts to disappear when NGC~6302 is observed in the Q-band. At longer wavelengths, where the spectral energy distribution (SED) peaks and the different solid state species are identified, the dust torus will be optically thin.

A simplified calculation verifies this: The dust is thought to be present in a torus, with a scale height $`h`$ from the equatorial plane and an outer radius $`R_{\mathrm{out}}`$. The inner radius is small compared to the outer radius and can be ignored. The torus is seen under an inclination angle of 45o (Bohigas 1994). The length of the line-of-sight through the disk is thus $`l=2h\sqrt{2}`$. Assuming that the density is constant throughout the disk, we find for the infrared optical depth along this line-of-sight:

$$\tau _{\mathrm{IR}}=Qn\sigma _{\mathrm{geom}}\mathrm{\hspace{0.17em}2}h\sqrt{2}$$
(8)
where

$$n=\frac{M_\mathrm{d}}{M_{\mathrm{grain}}}\frac{1}{2h\pi R_{\mathrm{out}}^{}{}_{}{}^{2}}$$
(9)
By adding all cold dust components listed in Table 2, the dust contained in the torus is determined to be $`M_\mathrm{d}=0.05`$ $`M_{}`$. The mass of a grain follows from the used grain size of 0.1 $`\mu `$m, and the average density of silicate grains. Using cgs units, we find that the torus is optically thin in the mid- and far-infrared in case

$$\frac{Q_{\mathrm{IR}}}{R_{\mathrm{out}}^{}{}_{}{}^{2}}<1.110^{36}\mathrm{cm}^2$$
(10)
For a typical value of $`Q_{\mathrm{IR}}`$ in the mid- and far-infrared of 0.0028 we find that for $`R_{\mathrm{out}}>5.010^{16}`$ cm, or $`3.410^3`$ AU, the dusty torus is optically thin at these wavelengths. This size is not unreasonable for a circumstellar torus, as it is comparable to the sizes of the dust shell shells around AGB stars, the PN progenitors. In addition, we can argue that the density distribution will not be flat, but probably behaves more like a power law. Then, most lines-of-sight through the torus will be optically thin even in case of a smaller torus. Therefore, we are confident that the torus is indeed optically thin in the mid- and far-infrared.

ii) We assumed that all grains have a size of 0.1 $`\mu `$m. In reality, the grains formed in the circumstellar environment will cover a wide range of sizes. However, we believe that we actually get a good estimate of the actual dust masses, because we use the mass absorption coefficients, and is thus independent from grain size. Using the mass absorption coefficients will lead to a direct determination of the dust mass, without calculating the mass of individual grains. The calcite and dolomite data presented in this work are mass absorption coefficients $`\kappa _\nu `$, the data used for diopside, enstatite and forsterite are measured in terms of absorption efficiency $`Q_\nu /a`$. For amorphous olivine, metallic iron and water ice we used optical constants, and derived absorption efficiency for grains with a size of 0.1 $`\mu `$m. Only the mass determination of these three species is probably somewhat affected by our assumption on the grain size.

iii) We assumed that all grains are spherical. Again an assumption that is probably not very realistic. The dominant dust component (amorphous olivine) must be in the form of non-spherical particles to explain the infrared-spectrum of AGB stars, and the same is true for metallic iron (Kemper et al. 2002a). The similarity in the spectral appearance of the solid state features detected in PNe and their direct progenitors, AGB stars, (e.g. Sylvester et al. 1999; Molster et al. 2002b), indicates that the dust grains in the two types of objects probably have similar properties. Therefore, we may conclude that the grains in NGC~6302 are very likely not homogeneous spheres, as these particles result in a shift in peak position with respect to non-spherical particles and inhomogeneous spheres. However, in our model the only purpose of the grain shape is to calculate the total mass of some of the dust components. The effect that grain shape has on light scattering is not taken into account, as our model uses mass absorption coefficients.

iv) We assumed that all grains are of homogeneous composition. Using mass absorption coefficients to calculate the emerging spectrum, we simply derived the total mass of a certain mineral found in the line of sight. It is impossible to discriminate between homogeneous grains, and grains consisting of two or more mineral species, therefore we cannot exclude the presence of grains of mixed composition. For instance, it is well known that volatiles, such as water vapour, condense in layered structures on grains (Jura & Morris 1985). We stress that the derived total masses are independent from the assumed homogeneity of grains.

v) We assumed that the density and temperature gradients are simple power laws of the distance to the central star. The complex structure of NGC~6302 indicates that this is probably far besides the truth, but since the optical depth is low, all infrared radiation emitted by the grains is received. Therefore not the density and temperature as a function of distance to the central star, but the mass-temperature relation of the different dust species becomes important, in a similar fashion as described by Bouwman et al. (2000). Our simulation of this relation by a power law density and temperature gradient seems to be in agreement with the observed spectrum.

vi) Finally we assumed that all dust species are found in the same temperature range. However, if the different dust components are not in thermal contact with each other, i.e. if they are present in separate – but co-spatial – grain populations, the difference in optical properties, notably in the UV and visual, causes a difference in the temperature profile (Bouwman et al. 2000; Kemper et al. 2001). The determination of the temperature range is based on the overall shape of the SED, and we forced the crystalline components – which are the carriers of the narrow features superposed on the broad continuum – to have the same temperature profile as the amorphous component. Since the dust emission strongly depends on the temperature of the grains, this will probably be our main source of error in our mass estimates. However, dropping the constraint that all dust species are found in the same temperature range, leads to a large increase in free parameters for our simple model.

Considering all the effects discussed here, we estimate that our absolute mass determinations are accurate within a factor of two, if the distance toward NGC~6302 is 910 pc. If we take into account the large uncertainty in the distance determinations, the error in the absolute dust mass determination further increases. The uncertainty in the various laboratory measurements introduces an additional uncertainty of a factor of two, leading to a total error of a factor of $``$3. This inaccuracy in the laboratory measurements can be inferred from comparison between published data. The crystalline silicate features published by Koike et al. (1999) are intrinsically much stronger (up to 5 times), than the same features in measurements by other groups (Servoin & Piriou 1973; Steyer 1974; Jäger et al. 1998). Consequently, the crystallinity, defined by the mass of enstatite and forsterite as a fraction of the total dust mass, lies within 5% – 25%, and the values given in Table 2 are thus just lower limits. The error in the relative dust mass fractions, i.e. the mass ratios between the identified minerals, is much better, as they all depend in the same way on our model assumptions. We therefore estimate the error in the relative dust mass fractions of the order of 50%. The error in the absolute and relative flux levels is small compared to the inaccuracies discussed here, and is therefore not taken into account.

In order to further approach reality, full radiative transfer calculations should be performed. Using spatial information derived from imaging and measurements of outflow velocities it is possible to reconstruct the geometry of the nebula. With 2-dimensional full radiative transfer calculations over a large wavelength range, from UV to far-infrared, it is possible to resolve the temperature and density distribution throughout the PN, for the different dust components. Also, the effect of grain sizes, grain shapes and mineralogical homogeneity then becomes apparent in the spectroscopy. This will lead to a much more accurate description of the properties of the dust in the nebula. However, we chose not to do these calculations for the following reasons: First, optical constants are not available for some of the identified dust components. Second, our knowledge about the geometry of the nebula is still too limited to justify the huge effort of 2-dimensional radiative transfer calculations. These calculations introduce new free parameters, and it will be very hard to constrain the model parameters in an unambiguous way with the limited information available. Spatially resolved spectroscopy of the mid- and far-infrared regions will certainly improve the image we have of NGC~6302. Finally, to obtain reliable mass estimates the distance toward NGC~6302 should be more accurately known.

### 4.2 The astronomical relevance of carbonates

In paper I and this work, we present the first extrasolar detection of carbonates. The detection is based on far-infrared features of carbonate at 92 $`\mu `$m and dolomite, which contributes to the 60 $`\mu `$m complex. Due to the low abundances, the strong resonances present at short wavelengths (see Fig. 6) do not stand out in the spectra, as many other species contribute significantly to the opacity at these wavelengths. This is not the case at the far-infrared wavelengths, therefore it is possible to detect even very small amounts of carbonates, beyond reasonable doubt.

There have been a number of searches for carbonate features in the mid-infrared. The previously unidentified infrared (UIR) feature at 11.3 $`\mu `$m seen in planetary nebula NGC~7027 was attributed to carbonates (Gillett et al. 1973; Bregman & Rank 1975). However, this identification became unlikely with the non-detection of the 6.8 $`\mu `$m carbonate feature in NGC~7027 (Russell et al. 1977b), and the lack of characteristic features due to carbonates in the 22–35 $`\mu `$m region of the spectrum (McCarthy et al. 1978). The UIR feature at 11.3 $`\mu `$m in HD~44179, also known as the Red Rectangle, was attributed to the carbonate magnesite (MgCO3) (Bregman 1977). However, Cohen et al. (1986) showed PAHs are the most likely carrier of the UIR feature at 11.3 $`\mu `$m seen in various astrophysical environments. Thus, the sole detection of the 11.3 $`\mu `$m feature is no longer considered evidence for the presence of carbonates.

Interesting environments to search for carbonates are star forming regions and young stars. Carbonates form easily on the surface of planets when liquid water is present, but other formation mechanisms cannot be excluded (Paper I). Russell et al. (1977a) have determined an upper limit to the carbonate/silicate mass ratio in the interstellar medium toward Orion. They find that the this mass ratio is at most 0.05, based on the detection limit of the 6.8 $`\mu `$m feature. The 6.8 $`\mu `$m feature is seen in absorption toward embedded protostars. The resemblance of the spectral appearance of the 6.8 $`\mu `$m feature in W33A and in Interplanetary Dust Particles (IDPs) has lead to the conclusion that they have the same carrier, most likely carbonates (Sandford & Walker 1985; Tomeoka & Buseck 1986), although Tielens et al. (1984) argue that it is due to hydrocarbons. Cox (1989) observes the 6.8 $`\mu `$m feature towards protostar AFGL~961 and claims that it is partially due to carbonates along with a contribution of hydrocarbons. Hints of the 11.3 and 13.4 $`\mu `$m carbonate features are also claimed to be present. The similarity between the 6.8 $`\mu `$m features in IDPs and in interstellar lines-of-sight is reported by Quirico et al. (2000) as well, presumably coming from the same carrier, although carbonates are rejected as a carrier. From the high resolution ISO SWS spectroscopy it becomes clear that the spectral shape of the 6.8 $`\mu `$m absorption feature observed towards protostars cannot be explained by carbonates, but a satisfactory alternative is still lacking (Keane et al. 2001). Recently, the $``$92 $`\mu `$m carbonate feature has been detected towards protostar NGC~1333-IRAS~4 (Ceccarelli et al., in prep.).

Although extraterrestrial carbonates are quite often found in meteorites and IDPs, they are actually not easy to detect by means of astronomical observations. The only claim still standing at the moment is the detection of features at 26.5 and 31 $`\mu `$m due to dolomite and calcite respectively in the ISO SWS spectrum of Mars (Lellouch et al. 2000). However, the far-infrared features of dolomite and calcite at $``$60 and $``$92 $`\mu `$m are not seen on Mars (Burgdorf et al. 2000).

As pointed out in Paper I, until now carbonates are believed to be formed through aqueous alteration. In this formation process, carbon dioxide (CO2) dissolves in liquid water and forms carbonate ions (CO32-). If cations like Ca2+ or Mg2+ are found in the solution as well, carbonates can be formed as a lake sediment when saturation is reached. The presence of carbonates is seen as evidence of planet formation, as an atmosphere containing carbon dioxide is required, as well as liquid water on the surface of the planet. The connection with planet formation is the justification for the search for carbonates in the circumstellar environment of young stars. Carbonates are ubiquitous in our own solar system, as is shown from the composition studies of meteorites and interplanetary dust particles. Carbonate-containing meteorites are considered tracers of the formation history of the solar system (e.g. Endress et al. 1996; Bischoff 1998). However, the detection of carbonates in PNe suggests that alternative formation mechanisms exist and that the presence of carbonates no longer provides direct evidence for planet formation (Paper I).

## 5 Summary

We have determined the composition and distribution of the dust in the planetary nebula NGC~6302. We found that a warm (100-118 K) and a cool (30-60 K) dust component are present. The cool dust component is located in a circumstellar torus of which the inner part effectively shields the UV/optical and near-infrared radition from the central star. The torus contains amorphous olivine, forsterite, clino-enstatite, water ice, diopside, dolomite and calcite. Outside the solar system, the carbonates dolomite and calcite are only seen in NGC~6302 and NGC~6537 (see also Paper I), and the formation mechanism of carbonates in these environments is not yet understood. Diopside is so far only found in NGC~6302 (Koike et al. 2000) and in two OH/IR stars (Demyk et al. 2000). Water ice, forsterite, enstatite and amorphous olivine are very common in the dust shells of evolved stars (e.g. Molster et al. 2002a, b, c). The circumstellar torus also contains the OH maser reported by Payne et al. (1988).

The optical depth in the polar regions is smaller, and therefore, these regions can be observed in the optical in scattered light from the central star. These regions contain the warm dust component, which thermally emits in the mid-infrared. The N-band images are dominated by thermal emission from the warm dust. From the ISO spectroscopy, constraints on the dust composition could be derived: predominantly metallic iron (or another featureless dust component), amorphous olivine, forsterite and clino-enstatite are detected. For water ice, diopside, calcite and dolomite we have derived upper limits for the mass fraction in the warm component. Only 0.3% of the total dust mass is contained in the warm component, which has a different mineralogical composition than the cool component. UV radiation can penetrate the polar regions and thus we are able to explain the observed PAH emission, which is also reported by Casassus et al. (2000).

In addition, the mass absorption coefficients of carbonate minerals are presented in this work. We studied calcite, dolomite, ankerite, aragonite siderite and magnesite.

###### Acknowledgements.

We thank M. Matsuura for her help in obtaining the TIMMI2 observations and A. de Koter for the useful discussions. We are grateful for the support by the staff of the ESO 3.6m telescope. LBFMW and FK acknowledge financial support from an NWO ’Pionier’ grant. This work was partly supported by NWO Spinoza grant 08-0 to E.P.J. van den Heuvel.


format: text

text:
# Ultrafast non-linear optical signal from a single quantum dot: exciton and biexciton effects

## 1 Introduction

Exciton dynamics experiments are attracting continuos interest because of their suitability to explore dephasing effects of single exciton and multiexciton complexes. Furthermore, recent proposals for solid state quantum computing systems, have stressed the importance of properly controlling multiexciton coherences. Recently, the non-classical behaviour of the light emitted from a SQD, photon antibunching in the fluorescence spectrum, has been observed in the ultrafast regime. Additionally, the search for single photon sources has triggered experimental interest in transient multiexciton coherences. All of these exciton based phenomena reflect the importance of interparticle interactions (electrons and holes) and their couplings to the environment, e.g, phonons. Both coherent and incoherent effects are experimentally detected. Incoherent effects are related to the system-bath interaction and from a theoretical point of view they are usually modeled within a Markov approximation. Coherent effects are associated with multiple particle correlations and can be observed in the subpicosecond time scale. Consequently, a detailed understanding of the coherence persistence at short times is of most importance in nanostructure systems.

Time resolved FWM experiments give signatures of particle-particle correlations and particle-environment interactions. In a typical FWM experiment, two simultaneous excitation pulses co- or cross-circularly polarized, with wave vectors $`\stackrel{}{k}_1`$ and $`\stackrel{}{k}_2`$, propagate onto the sample. A third laser pulse, with wave vector $`\stackrel{}{k}_3`$, is sent time delayed at $`t_2=T`$. For moderate excitation intensities, exciton-exciton interaction arises. The time evolution of the third order polarization $`P^{(3)}(t,T)=|P^{(3)}(t,T)|e^{i\mathrm{\Phi }(t,T)}`$ describes the properties of the diffracted light in the $`\stackrel{}{k}_3+\stackrel{}{k}_2\stackrel{}{k}_1`$ direction. At very short times after the excitation, particle correlations build up. These correlations have been observed using spectral interferometry techniques in bulk systems. It is also possible to study the polarization phase dynamics itself, $`\mathrm{\Phi }(t,T)`$. Moreover, other experiments for measuring $`P^{(3)}`$, like pump and probe schemes, have also been performed. In particular in Ref. a simple four level model has been used to get the non-linear optical response ($`\chi ^{(3)}`$) of a SQD. From the frequency spectrum of $`\chi ^{(3)}`$, the possibility of getting exciton entanglement in a SQD has been predicted. However, in that work exciton correlations were included within a Mean Field Approach (MFA).

Exciton dephasing times ($`T_2`$) have been well characterized in bulk and quantum well systems. In particular, ultrafast dynamics experiments, like time-resolved speckle analysis, have been used to get dephasing times information in solids. On the other hand, recent experiments in bulk systems show the build up of many particle correlations in the femtosecond scale of time. Recently, all of these dynamical phenomena have been observed for first time in Self Assembled QDs using transient FWM spectroscopy where the time integrated FWM signal oscillates with the biexciton binding energy. Besides that, dephasing mechanism from Optical Longitudinal (LO) phonons has been observed in CdSe Bulk and quantum dots systems, showing clearly non-markovian effects associated to the coupling between carriers and LO phonons. Due to its atomic like characteristics, a SQD offers the unique possibility of manipulating the number of particles and their Coulomb interactions (exciton-exciton (X-X)) with direct consequences on decoherence processes control. Thus, rather long dephasing times ($`T_2`$ 40 ps) have been reported.

A natural question arises when ultrafast scale of times are involved: what kind of quasiparticle (free electron-hole pairs, exciton or bound excitons) are dominant?. In order to answer this question we consider the Phase Space Filling (PSF) effect and X-X interactions on the same footing. Within MFA, the FWM signal shows contributions coming from both, PSF effects as well as renormalized X-X interactions. However, this approximation fails to explain correlations at very short times. Therefore, MFA should be improved by considering quantum fluctuation effects to explain transient FWM experiments. In particular, more elaborated theoretical approaches like the one reported in Ref., where a truncation in the hierarchy of polarization equations of motion to fifth order in the optical electric field is proposed, are shown to improve MFA results by including four particle correlations. Aditionally, some previous reported works show that in bulk systems a mean-field approach to calculate the FWM signal gives errors because the exact X-X correlations are neglected.

In the present work, we go beyond previous limitations like few level systems and MFA. We use a truncation scheme similar to one developed for bulk systems, in which the contribution from exciton and exciton-exciton effects can be split and treated on an equal footing to any desired order in the optical field. However, the truncation scheme fails to evaluate exactly the X-X Coulomb interaction due to the intensive computational work involved in its application to bulk systems. In order to get feasible results, this scheme has been used mapping the problem to a one-dimensional Hubbard model. In the present work we undertake this problem and we include exact X-X correlations in order to find the non-linear optical signal, intensity and phase dynamics in a realistic model of a SQD in the ultrafast regime. In Section 2, we briefly review the theoretical background on which our results are based. Section 3 describes the main results of this work. Our conclusions are summarized in Section 4.

## 2 Theoretical Model

Our starting point is the Hamiltonian for a system of electrons and holes in a two-dimensional parabolic quantum dot excited by external laser pulses

$`H=`$ $`{\displaystyle \underset{\nu ,s}{}}E_{\nu ,s}^ec_{\nu ,s}^{}c_{\nu ,s}+{\displaystyle \underset{\nu ,s}{}}E_{\nu ,s}^hh_{\nu ,s}^{}h_{\nu ,s}`$ (1)
$`+`$ $`{\displaystyle \underset{\nu ,\sigma (s,s^{})}{}}\mu _{\nu ,\sigma }^{eh}(t)c_{\nu ,s}^{}h_{\nu ,s^{}}^{}+{\displaystyle \underset{\nu ,\sigma (s,s^{})}{}}\mu _{\nu ,\sigma }^{he}(t)^{}h_{\nu ,s^{}}c_{\nu ,s}`$
$``$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{\nu _1,\nu _2,}{\nu _3,\nu _4}}{}}{\displaystyle \underset{s,s^{}}{}}\nu _1,\nu _2|V_{eh}|\nu _3,\nu _4c_{\nu _1,s}^{}h_{\nu _2,s^{}}^{}h_{\nu _3,s^{}}c_{\nu _4,s}`$
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{\nu _1,\nu _2,}{\nu _3,\nu _4}}{}}{\displaystyle \underset{s,s^{}}{}}\nu _1,\nu _2|V_{ee}|\nu _3,\nu _4c_{\nu _1,s}^{}c_{\nu _2,s^{}}^{}c_{\nu _3,s^{}}c_{\nu _4,s}`$
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{\nu _1,\nu _2,}{\nu _3,\nu _4}}{}}{\displaystyle \underset{s,s^{}}{}}\nu _1,\nu _2|V_{hh}|\nu _3,\nu _4h_{\nu _1,s}^{}h_{\nu _2,s^{}}^{}h_{\nu _3,s^{}}h_{\nu _4,s}`$
where $`c(h)_{\nu ,s}^{},c(h)_{\nu ,s}`$ creates and destroys one electron(hole) in the state $`\nu `$ (labeled by quantum numbers $`(n_{e(h)},m_{e(h)})`$) with spin $`s`$ ($``$ , spin up, $``$, spin-down). $`V_{ee},V_{eh}`$ and $`V_{hh}`$, denote electron-electron, electron-hole and hole-hole Coulomb interactions, respectively. Single particle energies for electrons (holes) are given by $`E_{\nu ,s}^{e(h)}=(n_{e(h)}+m_{e(h)}+1)\omega _{e(h)}`$, with confinement energies $`\omega _{e(h)}=1/m_{e(h)}l_{e(h)}^2`$ ($`\mathrm{}=1`$), where $`l_{e(h)}`$ is the parabolic confinement length size of the electron(hole). The set of laser pulses used to excite the SQD is described by the envelope amplitude $`(t)`$ and associated dipole moments $`\mu _{\nu ,\sigma }^{eh}`$ which take into account pulse polarizations. The spins of the electron ($`s`$) and hole ($`s^{}`$) determine the corresponding polarization index $`\sigma `$ i.e. $`\sigma =\sigma (s,s^{})`$. By diagonalizing the time-independent part of the Hamiltonian Eq.(1) ($`(t)=0`$), the energies and wave functions for one and two excitons in a SQD are obtained, from which the non-linear optical response is calculated.

Exciton and Biexciton wave functions are given respectively, in terms of non-interacting electron-hole pairs, as

$`|n_{s,s^{}}^X`$ $`=`$ $`{\displaystyle \underset{\nu _1,\nu _2}{}}\mathrm{\Psi }_{\nu _1,\nu _2}^n(c_{\nu _1,s}^{}h_{\nu _2,s^{}}^{})|0`$ (2)
$`|m^{XX}`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{\nu _1,\nu _2,}{\nu _3,\nu _4}}{}}\mathrm{\Psi }_{\nu _1,\nu _2,\nu _3,\nu _4}^m(c_{\nu _1,}^{}c_{\nu _2,}^{}\pm c_{\nu _1,}^{}c_{\nu _2,}^{})(h_{\nu _3,}^{}h_{\nu _4,}^{}\pm h_{\nu _3,}^{}h_{\nu _4,}^{})|0.`$
$`|0`$ is the SQD ground state and $`\mathrm{\Psi }_{\nu _1,\nu _2}^n`$ and $`\mathrm{\Psi }_{\nu _1,\nu _2,\nu _3,\nu _4}^m`$ are the exciton and biexciton amplitudes, respectively, which take properly into account the carriers spin configuration. In the following formulas we will stick to the notation of Eq.(2) where the index $`n`$ labels an exciton state, whereas the index $`m`$ is reserved for the labeling of biexciton states.

In discussing the time resolved spectra, we will be concerned with the situation in which the SQD is excited by very short laser pulses ($`\delta `$ functions) in resonance with the $`1s`$ heavy-exciton state. Therefore, a big number of single particle states has to be included in the calculation of one and two exciton spectra. The applicability of the $`\delta `$-pulses limit can be extended to typical accoustic periods when interaction with LO phonons can be neglected, compared with the deformation potential. However, it will be very important in II-VI semiconductors due to this high enough polar coupling. Our results are related with III-V systems, in which this coupling is low. In this way, when coherent phenomena coming from the oscillations of exciton or biexciton-phonon complexes are lost, incoherent phenomena aquires more importance and after this time it is well known that a phenomenological model of decoherence, as the one we will adopt here, can be used.

We perform our theoretical calculations following a truncation scheme  of the Hilbert space of one exciton and two interacting excitons. After a long algebra it is possible to give analytical expressions for the transient non-linear optical signal, which is proportional to $`P^{(3)}`$. It can be described as the sum of three terms: one representing the PSF ($`P^{PSF}`$) contribution, another one corresponding to the MFA term ($`P^{MF}`$) and a third one related to the X-X contribution ($`P^F)`$. This last term describes exact four-particle correlations and memory effects in terms of the so called force-force time correlation function. All these terms can be calculated for a SQD described by the Hamiltonian Eq.(1), giving place to expressions similar but not equal to those of Ref. (the differences are due to the confinement geometry of the dot, see the discussion in the next section). The resulting expressions for each one of these contributions for a SQD are given explicitly by

$$P_{n_1,\sigma _1}^{PSF}(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}=i\alpha _{n_1,\sigma _1}\underset{\genfrac{}{}{0pt}{}{n_0,n_2,n_3,}{\sigma _0,\sigma _2,\sigma _3}}{}\alpha _{n_0,\sigma _0}\alpha _{n_2,\sigma _2}^{}\alpha _{n_3,\sigma _3}^{}\mathrm{\Theta }(T)\mathrm{\Theta }(t)e^{(i\omega _{x,n_0}\mathrm{\Gamma })T}C_{n_0,\sigma _0;n_1,\sigma _1}^{n_2,\sigma _2;n_3,\sigma _3}$$
(3)
where $`\alpha _{n,\sigma }`$ is related to the exciton wave function at zero relative distance, $`C_{n_0,\sigma _0;n_1,\sigma _1}^{n_2,\sigma _2;n_3,\sigma _3}`$ is the PSF parameter (for notation details, see Ref.) and $`\mathrm{\Theta }(x)`$ denotes the step function ($`\mathrm{\Theta }(0)`$=1). The MFA contribution can be written as

$`P_{n_1,\sigma _1}^{MF}(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}`$ $`=`$ $`i\alpha _{n_1,\sigma _1}{\displaystyle \underset{\genfrac{}{}{0pt}{}{n_0,n_2,n_3,}{\sigma _0,\sigma _2,\sigma _3}}{}}{\displaystyle \frac{\alpha _{n_0,\sigma _0}\alpha _{n2,\sigma _2}^{}\alpha _{n3,\sigma _3}^{}e^{(i\omega _{x,n_0}\mathrm{\Gamma })T}\beta _{n_0,\sigma _0;n_1,\sigma _1}^{n_2,\sigma _2;n_3,\sigma _3}}{\omega _{x,n_0}+\omega _{x,n_1}\omega _{x,n_2}\omega _{x,n_3}+2i\mathrm{\Gamma }}}\times `$ (4)
$`\times \left[\mathrm{\Theta }(T)\mathrm{\Theta }(t+T)A_x(t,T)+\mathrm{\Theta }(T)\mathrm{\Theta }(t)A_x(t,0)\right]`$
where

$$\beta _{n_0,\sigma _0;n_1,\sigma _1}^{n_2,\sigma _2;n_3,\sigma _3}=\underset{m}{}(\omega _{xx,m}\omega _{x,n_0}\omega _{x,n_1})B_{n_0,n_1}^{\sigma _0,\sigma _1}(m)B_{n_2,n_3}^{\sigma _2,\sigma _3}(m)$$
(5)
and

$$A_x(t,T)=e^{(2\mathrm{\Gamma }+i(\omega _{x,n_0}+\omega _{x,n_1}\omega _{x,n_2}\omega _{x,n_3}))t}e^{(2\mathrm{\Gamma }i(\omega _{x,n_0}+\omega _{x,n_1}\omega _{x,n_2}\omega _{x,n_3}))T}.$$
(6)
Similarly, the X-X contribution is given by

$`P_{n_1,\sigma _1}^F(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}=i\alpha _{n_1,\sigma _1}\left({\displaystyle \underset{m}{}}P_{n_1}(m)\right)P_{n_1,\sigma _1}^{MF}(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}`$ (7)
where

$`P_{n_1}(m)={\displaystyle \underset{\genfrac{}{}{0pt}{}{n_0,n_2,n_3,}{\sigma _0,\sigma _2,\sigma _3}}{}}\alpha _{n_0,\sigma _0}\alpha _{n_2,\sigma _2}^{}\alpha _{n_3,\sigma _3}^{}B_{n_0,n_1}^{\sigma _0,\sigma _1}(m)B_{n_2,n_3}^{\sigma _2,\sigma _3}(m){\displaystyle \frac{e^{(i\omega _{x,n_0}\mathrm{\Gamma })T}(\omega _{xx,m}\omega _{x,n_0}\omega _{x,n_1})}{\omega _{xx,m}\omega _{x,n_0}\omega _{x,n_1}i\mathrm{\Gamma }_{xx}}}\times `$
$`\times [(\mathrm{\Theta }(T)\mathrm{\Theta }(t+T)A_{xx}(t,T)+\mathrm{\Theta }(T)\mathrm{\Theta }(t)A_{xx}(t,0)]`$ (8)
with

$$A_{xx}(t,T)=e^{(\mathrm{\Gamma }_{xx}+i(\omega _{x,n_0}+\omega _{x,n_1}\omega _{xx,m}))t}e^{(\mathrm{\Gamma }_{xx}i(\omega _{x,n_0}+\omega _{x,n_1}\omega _{xx,m}))T}.$$
(9)
In Eqs.(5) and (7), $`m`$ runs over all (bound and unbound) biexciton states. Bound biexcitonic states are those for which its energy lies below two times the ground state energy of an exciton. Aditionally, $`\beta _{n_0,\sigma _0;n_1,\sigma _1}^{n_2,\sigma _2;n_3,\sigma _3}`$ and $`P_{n_1,\sigma _1}(m)`$ take into account the exciton and biexciton transition weights on the signal produced by the $`n_1`$ exciton state. The exciton transition weights $`B_{n_i,n_j}^{\sigma _i,\sigma _j}(m)`$ are defined to be the matrix elements $`0|\widehat{B}_{n_i,\sigma _i}\widehat{B}_{n_j,\sigma _j}|m^{XX}`$ of the product of two $`n^{th}`$ exciton destruction operators $`\widehat{B}_{n,\sigma }`$ between the exciton ground state $`0|`$ and a biexciton state $`|m^{XX}`$. Finally, all we need is: exciton, $`\omega _{x,n}`$, and biexciton, $`\omega _{xx,m}`$, energies and their respective wave functions, which can be numerically obtained to any desired precision in the SQD case. From these ingredients, the exciton transition weights can be calculated as $`B_{n_0,n_1}(m)=0|\widehat{B}_{n_0}\widehat{B}_{n_1}|m^{XX}`$ where the exciton destruction operator is $`\widehat{B}_{n,\sigma }=_{\nu _1,\nu _2}\mathrm{\Psi }_{\nu _1,\nu _2}^nc_{\nu _1,s}h_{\nu _2,s^{}}`$ ($`\sigma =\sigma (s,s^{})`$) and $`\mathrm{\Psi }_{\nu _1,\nu _2}^n`$ is the $`n`$-th exciton amplitude (cf. Eq.(2)). One of the main results of our paper is that, in contrast to an extended system, for QDs the term $`P_{n_1,\sigma _1}^{MF}`$ does not vanish, *even* in the case of cross-polarized (CP) excitation. If the dot confinement energies are reduced, in the CP case $`P_{n_1,\sigma _1}^{MF}`$ decreases in magnitude, approaching zero in the limit of zero confinement (extended system). This can be seen by direct evaluation of the parameter $`\beta `$ in Eq.(5), that gives information on the strength of the signal as predicted by MFA. Although, the $`\beta `$ term contains information about the renormalized biexciton energy, this term is not responsible for oscillations and only contributes with an spectral weight given by the $`B_{n_i,n_j}^{\sigma _i,\sigma _j}`$ changing the intensity of the FWM signal.

The total dynamics of excitons and biexcitons comprises the MF and the exact X-X contributions. The temporal dynamics dependence of $`P^{MF}`$ (apart from an exponential damping) is contained in Eq.(6) where only contributions from differences of exciton energy levels appear, which are negligible in the case of resonant excitation. On the other hand, the temporal dependence of the X-X contribution is dominated by energy differences between two exciton states and one biexciton state, e.g. Eq.(9). In particular, in the CP case, a bound biexcitonic state is present, causing oscillations in the signal with a frequency corresponding to the energy of the bound state.

Exciton and biexciton dephasing rates associated to non-radiative mechanisms are described (phenomenologically) by $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }_{xx}`$, respectively. In the following we will concentrate only in the lowest exciton state, i.e. $`n_1=n_2=n_3=n_4=1s`$ with $`\mathrm{\Gamma }_{xx}=2\mathrm{\Gamma }`$. In our model, the dephasing time $`T_2`$ is defined as the inverse of $`\mathrm{\Gamma }_{xx}`$.

## 3 Results and discussion

In order to calculate the FWM signal for a typical self-assambled SQD, we consider a $`In_{0.5}Ga_{0.5}As`$ quantum dot with $`ϵ=12.5`$, $`m_e^{}=0.068`$, $`m_h^{}=0.2`$, $`\omega _e=20`$ meV and $`\omega _h=5`$ meV. To guarantee a good convergence for the spectrum of zero total angular momentum excitons and biexcitons, we diagonalize a 72x72 Hamiltonian matrix for excitons and a 3250x3250 Hamiltonian matrix for biexcitons. The exciton spin is determined by the circular polarization of the excitation pulses. Further in this work we discuss only the cross polarization excitation case. With this set of parameters the biexciton binding energy is -1.27 meV

As discussed in Section 2, the FWM signal clearly depends on $`T`$, on the exciton and biexciton energies as well as on the transition spectral weights. Excitation pulses create coherent excitons whose dephasing time ($`T_2`$) can be monitored by the delayed pulse. As the time delay increases, the FWM signal decreases from which $`T_2`$ can be obtained. With these ingredients, and considering that more than one electronic excitation can be induced by the laser and four particle correlations occur for times lesser than $`T_2/4`$, we expect that interference phenomena appear in the ultrafast time scale. Therefore, the coherent FWM signal and the phase dynamics give clear signatures of exciton-exciton interactions in the ultrafast regime. In the following we will not be concerned with $`T_2`$ determination, which is to be taken as a free parameter. Pulses with zero time delay, $`T=0`$, will be considered.

In the weak intensity regime both PSF and X-X effects determine the nonlinear optical properties. We will concentrate on these two contributions. PSF contributions to the FWM signal are only important in the parallel spin exciton case, where X-X effects are negligible, because a bound biexciton state does not exist at all. By contrast, for the cross polarization situation we are interested in, PSF effects are vanishing while X-X effects become dominant. It is important to remark that in extended and isotropic systems the MFA gives a zero contribution to the FWM signal, because only excited excitons with zero center of mass momentum contribute, i.e.,($`\stackrel{}{q}=0`$). However, in SQD this condition can be relaxed due to the fact that the parabolic potential absorbs the incident momentum and excitons with center of mass momentum different from zero can be excited. This point becomes particularly clear by considering the mean field parameter $`\beta `$ defined in Eq.(5). Whereas for extended systems this parameter is exactly zero , for a SQD it becomes a function of the confinement energies $`\omega _{e(h)}`$. The dependence of $`\beta `$ on $`\omega _e`$ is depicted in Fig. (1). We consider a constant relationship between electron and hole confinement lengths, i.e, $`l_e=l_h`$. It can be seen that in the limit of small confinement, $`i.e.`$ pure two-dimensional systems ($`\omega _{e,h}`$=0), $`\beta `$ tends to zero, as it should be. This means that the usual condition of ($`\stackrel{}{q}0`$) for the carriers Coulomb interaction in extended systems is broken in confined systems. We have tested our results taking a huge set of excitonic states and we found that higher energy states give a negligible contribution to $`\beta `$.

The phase-dynamics and the non-linear optical signal are calculated for different decay rates. In order to better clarify, why MFA is not adequate to describe non-linear optical signals at small times, we plot in Fig. 2, separately, the amount of particle correlation contribution to the FWM signal as described within MFA (Eq.4) and exact X-X effects (Eq.7), which include memory (non-Markovian) effects, for $`\mathrm{\Gamma }_{xx}=0.125`$ meV. The usual and simplest approximation for this last term yields to the Markov approximation and MFA results. By contrast, in this work we have performed a full numerical evaluation of X-X effects by summing over all biexciton states. Hence, we plot the absolute values of both, the real part of $`P_{n_1,\sigma _1}^{MF}(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}`$ (MFA term) and the real part of $`P_{n_1,\sigma _1}^F(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}`$ (X-X term). The MFA result grows monotonically with $`t`$, whereas the X-X result shows beating oscillations with a period corresponding to $`T_{xx}=2\pi \omega _{xx,opp}^15ps`$, the biexciton binding energy. As can be seen, non-Markovian effects arising from exciton-exciton correlations are the main source of discrepancies between MFA results and the exact ones. This means that ultrafast spectroscopy in SQDs should register memory effects associated to four particle correlations.

In the inset of Fig. 2, we show the spectral weight $`B(m)`$ as a function of biexciton energy levels (the zero of energy corresponds to twice the single $`1s`$-heavy-exciton energy). We have only plotted those levels that have total angular momentum $`L_z=0`$. The bound biexciton level is clearly seen at negative energies. Although, unbound biexciton states have non-negligible spectral weights the bound biexciton still has a comparable weight. As can be seen in Fig. 2 MFA fails to describe properly four particle interactions at very short times, where particle-particle correlations are important. This term does not oscillates and only average the exact result. Besides that, the contribution in Eq.7, takes into account a higher degree of correlation going beyond MF and gives oscillations at very short times. This oscillations provide information about strong correlations between excitons. Thus, the bound biexciton dynamics is directly linked to non-Markovian effects. By contrast, in the long time limit ($`t>T_2`$), memory effects are lost and MFA results approach to the real part of the first term in Eq.7. The FWM signal becomes dominated by non-radiative dephasing effects, producing in this way a simple exponential decay.

In order to assess how the main features described above are reflected on the optical dynamical characteristics, we plot the non-linear FWM intensity $`|P^{(3)}(t,0)|`$ in Figs. 3(a) and 3(c), and the dynamical phase $`\mathrm{\Phi }(t,0)`$ in Figs. 3(b) and 3(d), for dephasing rates, $`\mathrm{\Gamma }_{xx}=T_2^1=`$0.125 meV and 0.5 meV, respectively. To better illustrate the results, a comparison is made between MFA and X-X results on the same plots. For a low dephasing rate ($`T_2`$ 8 ps), strong oscillations are clearly seen in Figs. 3(a) and (b) for $`t<T_2`$. The X-X interaction, induces a $`\pi `$ phase shift with respect to the incident electric field. The phase of the $`P^{(3)}`$ signal oscillates during a typical time $`T_2`$. FWM signal shows a strong peak at $`tT_2/4`$ whereas it should have a maximum at $`tT_2`$ from MFA. For high dephasing rates ($`T_2`$ 2 ps), Figs. 3(c) and (d), the attenuation of both the intensity and phase of the third order optical signal is evident. At this time scale, the $`T_2`$ dephasing time is shorter than the time for which coherent effects could be seen. It is worthnoting that while the phase shows a similar starting behaviour as compared with the low dephasing case, the FWM signal intensity is drastically reduced by roughly a factor of five (notice the change of vertical scale).

Most importantly, the phase dynamics exhibits novel features at small times when exciton-exciton correlations are considered. By contrast, MFA results yield to a real $`P^{(3)}`$ for any time. Therefore, no phase dynamics is observed within MFA. In particular, X-X effects present a correlation between the first FWM signal minimum and a $`\pi `$ jump in its phase. For long times the phase goes to a constant value. However, after the extinction of the FWM intensity, there is no longer a clear physical meaning for this phase. Hence, oscillations in the phase dynamics, with a beating frequency controlled by the biexciton binding energy, could bring important information about exciton correlations in SQDs. In contrast with higher dimensional systems (quantum wells and bulk samples), a SQD provides us with an adequate system where exciton and biexciton discrete energy spectra can be tailored by changing the confinement potential. In this way, optimal conditions to enhance phase memory could be explored.

## 4 Conclusions

In summary, we have studied the non-linear optical response from multiexciton complexes in a SQD. For cross-polarized excitation, the creation of bound biexcitons is possible. They dominate the non-linear polarization dynamics in the low density regime. The phase dynamics and the FWM signal in the ultrafast scale of time have been obtained. By including all the numerically determined biexciton states, exact exciton-exciton correlations are evaluated allowing to go beyond simple MFA results. In particular we found that the phase dynamics and the FWM intensity oscillate with a period which coincides with the inverse of bound biexciton frequency. We demonstrated clearly that the usual Mean Field theory, which is closely related to the Markovian approximation, does not describe correctly exciton correlations at very small times. However, our results show that at long times, MFA results for both intensity and phase dynamics, are close to those produced by a more realistic calculation.

Acknowledgents FJR and LQ acknowledge partial support from COLCIENCIAS (Colombia) project No.1204-05-10326/11408, Banco de la República (Colombia) and from the EPSRC-DTI LINK project (UK).

## References

Figure Captions

Fig. 1: MFA contribution ($`\beta `$ term) as a function of the electron and hole confinement frequencies, for $`\mathrm{\Gamma }_{xx}`$=0.125 meV.

Fig. 2: Absolute values of the MFA contribution (real part of $`P_{n_1}^{MF}(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}`$) and X-X contribution ($`P_{n_1}^F(t,T)e^{(i\omega _{x,n_1}+\mathrm{\Gamma })t}`$) to the FWM signal, for $`\mathrm{\Gamma }_{xx}`$=0.125 meV. Inset: Spectral weight of the biexciton wave function. Bound biexciton energy=-1.27 meV. The long-dashed line denotes the zero energy.

Fig. 3: $`\mathrm{\Gamma }_{xx}`$=0.125 meV: (a) Time resolved FWM signal; (b) Phase dynamics. $`\mathrm{\Gamma }_{xx}`$=0.5 meV: (c) Time resolved FWM signal; (d) Phase dynamics.