text:
# The mineralogy, geometry and mass-loss history of IRAS 16342-3814 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) with the participation of ISAS and NASA. The SWS is a joint project of SRON and MPE. Also based on observations obtained at the European Southern Observatory (ESO).

## 1 Introduction

As a Main Sequence star with an initial mass between $`1`$ and about $`8\mathrm{M}_{}`$ ages, it eventually evolves into an Asymptotic Giant Branch (AGB) star. An AGB star is a luminous cool giant star in which a dust driven wind expels matter from its envelope at a high rate ($`10^7`$ to $`10^4\mathrm{M}_{}\mathrm{yr}^1`$). Mass-loss increases over time as the AGB star evolves. This eventually ends in an episode of extremely high mass-loss, the so-called superwind (SW) phase, in which the mass-loss rate, $`\dot{M}`$, exceeds $`10^5\mathrm{M}_{}\mathrm{yr}^1`$. The SW rapidly strips the star of virtually its entire envelope. The high mass-loss during the AGB phase (the total AGB phase lasts some $`10^6\mathrm{yr}`$) creates a circumstellar envelope which may completely obscure the star. When the SW phase stops, the central star starts to shrink in radius at constant luminosity, getting hotter over time. The star is now said to be in the post-AGB phase. The old circumstellar envelope, which initially obscures the central star, becomes optically thin at optical wavelengths after $`1001000\mathrm{yr}`$ because it is still expanding. Depending on the rate at which the central star evolves, the star may start to photoionize the material that was ejected during the AGB phase, creating a planetary nebula (PN). In this case the phase between the end of the AGB and the birth of the PN is called the proto-planetary nebula (proto-PN) phase.

In this paper we study the proto-PN IRAS 16342-3814 (hereafter: IRAS 16342; also OH 344.1+5.8 or ”the water-fountain nebula”). The evolutionary status of IRAS 16342 as a proto-PN is examined in detail by Sahai et al. (1999, hereafter: SAH99). They report the presence of an asymmetrical bipolar reflection nebula in Hubble Space Telescope (HST) images at $`0.55`$ and $`0.80\mu `$m. The two lobes of the nebula, where the eastern lobe is significantly redder than the western lobe, are separated by a dark equatorial waist. The nebula is illuminated by starlight escaping through the polar holes in this ‘optically thick, dense, flattened cocoon of dust’. The central star is completely obscured by the dust cocoon (hereafter called ‘the torus’). The individual symmetry axes of the lobes do not coincide, although their orientation and location with respect to the dark waist are centrally point symmetric. OH maser line data at 1612, 1665 and 1667 MHz overlaid on the $`0.80\mu `$m image suggest that the nebular axis is tilted such that the western lobe is nearest to us. Faint halos are seen around the lobes. The halo regions are illuminated by scattered light from the lobes, and possibly by heavily reddened, direct starlight (SAH99).

IRAS 16342 has been observed as part of a mid-infrared (Meixner et al. 1999) and optical imaging survey (Ueta et al. 2000). In the mid-infrared survey IRAS 16342 was observed at $`9.8`$, $`12.5`$, $`17.8`$ and $`20.6\mu `$m. At these wavelengths the star was classified as a core/elliptical source, i.e. a source with an unresolved, very bright core that is surrounded by a lower surface brightness elliptical nebula. At all mid-infrared wavelengths the source was found to extend northeast-southeast with the same position angle for all wavelengths.

Early studies on IRAS 16342 revealed the presence of extremely high velocity outflows in both $`\mathrm{H}_2\mathrm{O}`$ ($`130\mathrm{km}\mathrm{s}^1`$) and OH ($`70\mathrm{km}\mathrm{s}^1`$) maser line emission (Zuckerman & Lo 1987; Likkel & Morris 1988). Likkel & Morris interpreted the $`\mathrm{H}_2\mathrm{O}`$ masers as streams or clumps of molecular gas forced out at the polar axes, and believed the OH masers to originate from lower velocity material at intermediate latitudes. Likkel et al. (1992) suggested a binary star model for IRAS 16342, as this could account for both the bipolarity and the high outflow velocities of the masers.

Radiative transfer calculations on IRAS 16342 have been performed by Guertler et al. (1996) and Zijlstra et al. (2001). The latter authors found an inner radius of $`1^{\prime \prime }`$ and an outer radius of $`1.5^{\prime \prime }`$ for the torus, in agreement with the HST observations. The mass of the torus was found to be $`0.1\mathrm{M}_{}`$. The radiative transfer calculations by both groups suggest the presence of a (very cool) outer shell.

In this paper we present the $`2200\mu `$m Infrared Space Observatory (ISO) spectrum and infrared ISAAC and TIMMI2 images ($`3.8`$ to $`20\mu `$m) of IRAS 16342. The aim of this paper is to study the mineralogy, morphology and mass-loss history of IRAS 16342. The paper should also add to our understanding of proto-PNe in general and to the injection of matter into the ISM by AGB stars. The latter is important since AGB stars are among the main ‘factories’ of interstellar dust in galaxies. In Sect. 2 the data and data reduction procedures are described. In Sect. 3 we present and discuss the ISO spectrum. Based on the spectrum, we study the mineralogy and mass-loss history of IRAS 16342. The ISAAC and TIMMI2 images of IRAS16342 are presented in Sect. 4 and will be used to study the morphology of IRAS 16342 and set limits on the grain sizes in the torus and the lobes. The conclusions are listed in Sect. 5.

## 2 Data and data reduction

### 2.1 ISO data

The $`2200\mu `$m spectrum was obtained with ISO (Kessler et al. 1996), using the Short Wavelength Spectrometer (SWS; de Graauw et al. 1996) for the $`245\mu `$m region and the Long Wavelength Spectrometer (LWS; Clegg et al. 1996) for the $`45200\mu `$m region. The two spectra were reduced independently. For both the SWS and LWS reduction, we used pipeline version OLP 10 as the starting point. In both cases, individual subbands were cleaned from glitches, flatfielded (to the mean of all detectors), sigmaclipped (using the default values $`\sigma =3`$ for SWS and $`\sigma =2.5`$ for LWS) and rebinned ($`\lambda /\mathrm{\Delta }\lambda =500`$ for SWS and $`\mathrm{\Delta }\lambda =0.068\mu `$m for LWS). The subbands were scaled (to band 1 and/or 2 for SWS and detector 5 for LWS) to obtain a continuous spectrum. Finally, the LWS spectrum was scaled to the SWS spectrum to get one continuous spectrum. Without scaling, the flux levels of the two spectrographs matched at $`45\mu `$m to within $`12\%`$. For LWS there was an off source spectrum available to correct for interstellar flux contributions. The flux levels in this spectrum were low ($`10\mathrm{Jy}`$), showing that the on-source spectrum has no significant background contribution. We have compared our reduction with the available IRAS data. Between $`14`$ and $`25\mu `$m, the ISO spectrum overlaps the IRAS Low Resolution Spectrum (LRS) to within $`5\%`$, while between $`12`$ and $`14\mu `$m the IRAS LRS flux is about $`25\%`$ higher than the ISO flux. The IRAS Point Source Catalog (PSC) fluxes match those derived from the ISO data to within $`26\%`$ (see Table 1).

### 2.2 ISAAC data

The IRAS 16342 $`3.8\mu `$m (L-band) image was obtained with the Infrared Spectrometer And Array Camera (ISAAC) 111http://www.eso.org/instruments/isaac/index.html at the Very Large Telescope (VLT) at Paranal, Chile, on August 6, 2000, using the Santa Barbara Research Center Aladdin 1024$`\times `$1024 pixel InSb array. The pixel scale of this detector is 0.071 ′′ per pixel. The total field of view is 72′′ $`\times `$ 72′′. The L-band filter has a central wavelength of 3.78 $`\mu `$m and a width of 0.58 $`\mu `$m. The optical seeing during the observations varied between 1.2′′ and 2.0′′.

The image was obtained using a 10′′ chop throw (perpendicular to the long axis of the source), while nodding was performed. In addition to the chopping, the source position was also jittered. The total integration time was 5.0 minutes, not considering the negative image in the chopped image. Only the positive image was used, since the PSF of the negative image was systematically wider. We used two photometric calibration stars: Y3501 and HR6736. The L-band magnitudes of these stars were measured by MSSSO photometry (McGregor 1994) and van der Bliek et al. (1996), respectively. As a reference for the point spread function (PSF), a nearby star, USNO 450$``$24728679 ($`m_\mathrm{L}=7.0\mathrm{magn}`$) was used. The obtained PSF was used to deconvolve the ISAAC image of IRAS 16342. The FWHM of the PSF was calculated to be 0.35′′. The reduction was done using eclipse package version 4.1.2 (Devillard 1997), and DAOPHOTO under IDL.

### 2.3 TIMMI2 data

IRAS 16342 images at $`7.9`$ (narrow band), $`9.8`$ (narrow band), $`10.6`$ (broad band) and $`20.0\mu `$m (broad band) were observed with the Thermal Infrared Multi Mode Instrument (TIMMI2) at the ESO 3.60m telescope at La Silla, Chile, on August 11, 2001, using a Raytheon 240$`\times `$320 pixel AsSi array. The pixel scale of this detector is 0.202′′ per pixel. The recorded optical seeing during the observations varied between 0.92′′ and 2.0′′.

The images were obtained using a 15′′ chop throw in the north-south direction while nodding was performed. As a reference for the PSF, HD 124897 ($`=\alpha `$ Boo) and HD 178345 were used. We combined individual chopped images into one image to improve the signal-to-noise. This was done by shifting the positive and negative images in each chopped frame to a prescribed position and then adding or subtracting them. Frames of poor quality were omitted. Finally, we deconvolved the resulting image, using either HD 124897 or HD 178345 as PSF. Due to non-photometric weather during the observations, the absolute flux calibration of the TIMMI2 images proved difficult and the results uncertain. Therefore, further attempts to do an absolute flux calibration on the images were not made. Moreover, we are mainly interested in the morphology of IRAS 16342 at different wavelengths, which does not require an accurate (absolute) calibration of the flux. We will use relative fluxes in the individual images only.

## 3 The ISO spectrum

### 3.1 Description of the spectrum

In the top panel of Fig. 1 we present the $`2200\mu `$m ISO spectrum of IRAS 16342 (solid line). A possible continuum for the spectrum, which we will discuss below, is also given (smooth solid line). The second panel shows the spectrum of the OH/IR star AFGL 5379, studied in detail by Sylvester et al. (1999). The continuum for AFGL 5379 is taken from that paper. The $`10`$ and $`20\mu `$m amorphous silicate features are clearly visible. These absorption features are observed in many AGB stars, their strength depending on the optical depth of the dust, which is directly related to the mass-loss rate. The absorption features become stronger, and the spectrum redder, as $`\dot{M}`$ increases (e.g. Bedijn, 1987). IRAS 16342 also shows the amorphous silicate features. The one at $`20\mu `$m was already detected in the IRAS LRS spectrum and fitted by Zijlstra et al. (2001). Compared to AFGL 5379, the amorphous silicate features in the spectrum of IRAS 16342 might seem to be absent at first glance! The spectrum of IRAS 16342-3814 is that of an extremely reddened OH/IR star. The third panel in Fig. 1 shows the spectrum of HD 161796. HD 161796 is another post-AGB object which has been studied in detail by Hoogzaad et al. (2002) and Meixner et al. (2002). The continuum for HD 161796 is taken from Hoogzaad et al. (2002).

A first inspection of the spectrum of IRAS 16342 suggests the presence of the characteristic features of crystalline silicates and crystalline water ice, though up to $`45\mu `$m they appear in absorption instead of in emission. To illustrate this, the continuum shown in the top panel of Fig. 1 has been drawn accordingly (smooth solid line). Following Molster et al. (2002b), we represented the continuum by a spline-fit and required it to be smooth and to maximize the continuum flux in both $`F_\nu `$ and $`F_\lambda `$. We emphasize that this continuum is used only to enhance the visibility of sharp features in the spectrum and has no physical meaning. In the lower panel of Fig. 1 we plot the continuum subtracted spectra of IRAS 16342 (solid line) and HD 161796 (dashed line). The continuum subtracted spectrum of HD 161796 is dominated by crystalline silicate and crystalline water ice emission features (Hoogzaad et al. 2002). A comparison of the suggested absorption features in IRAS 16342 with the emission features in HD 161796 shows an impressive agreement, indeed identifying the IRAS 16342 features as due to crystalline silicates and water ice, and justifying our original choice of the continuum. At wavelengths longwards of approximately $`45\mu `$m the features in the spectrum of IRAS 16342 are again in emission. Close to the peak of the SED, where the transition between absorption and emission occurs, the shape of the spectral features is less certain.

### 3.2 The mineralogy of IRAS 16342

Table 2 lists the wavelength positions and identifications of the features found in the spectrum of IRAS 16342. We identify amorphous silicates, forsterite ($`\mathrm{Mg}_2\mathrm{SiO}_4`$), diopside ($`\mathrm{MgCaSi}_2\mathrm{O}_6`$), crystalline water ice, and possibly clino-enstatite ($`\mathrm{MgSiO}_3`$) as important solid state components of the circumstellar matter around this object. We also find \[OI\] and \[CII\] lines, indicating low density gas. The off source spectrum suggests that the \[CII\] line is largely interstellar, although a circumstellar contribution (up to $`25\%`$ in peak strength) may still be present. The \[OI\] line seems to be interstellar. After subtraction of the off source spectrum, in which the \[OI\] line is clearly seen, there is no clear evidence for a residual circumstellar component.

Fig. 2 shows part of the continuum subtracted spectrum of IRAS 16342 in more detail, together with the mass absorption coefficients for forsterite, clino-enstatite and diopside (Koike et al. 1999). Also shown is the absorbance of crystalline water ice measured at 50 K (Smith et al. 1994). The wavelength positions of the identified forsterite, diopside and crystalline water ice features are indicated with short vertical lines. We also indicate the positions of clino-enstatite bands. Due to blending with forsterite, and given that the bands are fairly weak, we cannot be certain that clino-enstatite is present. However, in most evolved stars both forsterite and enstatite are found together (Molster et al. 2002a), and we suggest this is also the case in IRAS 16342. The $`18.5\mu `$m amorphous silicate feature, the \[OI\] line and the unidentified $`48\mu `$m feature (e.g. Molster et al. 2002a) are also indicated.

The solid state materials found in the circumstellar environment of IRAS 16342 are also found around other (high mass-loss), oxygen rich, evolved stars (e.g. Koike et al. 2000; Molster et al. 2002a). The very red ISO spectrum suggests that the bulk of the material around IRAS 16342 must be cold. A low temperature may also be implied by the $`69\mu `$m feature of forsterite. This intrinsically weak feature is detected either when the forsterite is cold or very abundant. The feature at $`48\mu `$m has been found in the spectra of many sources (e.g. Molster et al. 2002a). Its nature is unclear at this time. FeSi has been proposed as a possible carrier (Ferrarotti et al. 2000), but this identification remains debatable (Molster et al. 2002a). The $`48\mu `$m feature is seen in emission, superposed on the $`45\mu `$m feature of crystalline water ice, which is in absorption. If real, this suggests that the carrier of the feature has a different spatial distribution than the crystalline water ice, otherwise both materials would be seen in absorption. An alternative explanation is that the feature is an artifact in the responsivity function of the ISO LWS detector. More study will be needed to explore this issue.

### 3.3 The Mass-loss history of IRAS 16342

No other OH/IR star known to us displays a spectrum reddened to the extent as is seen in IRAS 16342. Also, no other OH/IR star is known to have crystalline silicate features in absorption at wavelengths up to $`45\mu `$m. Though Sylvester et al. (1999) showed that crystalline silicate absorption features are present in OH/IR stars, they did not detect them longwards of $`25\mu `$m. The behaviour of the spectrum and the crystalline features in it is the result of a very high optical depth in the dust shell, which suggests that the star must have an extremely large mass-loss rate.

We modeled the spectrum of IRAS 16342 using the radiative transfer programme modust (Bouwman 2001). We did not try to make a detailed fit of the spectrum, but rather wanted to get an order of magnitude estimate on the mass-loss rate. We assumed an M9III AGB central star ($`T_{\mathrm{eff}}\mathrm{\hspace{0.17em}2670}\mathrm{K}`$ and $`R_{}=372\mathrm{R}_{}`$) surrounded by a spherical circumstellar envelope (composed of olivine and iron particles, following Kemper et al. (2002)) created by mass-loss from the central star. We derived $`\dot{\mathrm{M}}10^3\mathrm{M}_{}\mathrm{yr}^1`$. In this case the best match between the peak position of the synthetic spectrum and that of the data was obtained. Also, the lack of near infrared ($`210\mu `$m) flux was best modeled in this case.

The crystalline silicate absorption features in the ISO spectrum require that the envelope is optically thick up to about $`45\mu `$m. This implies that the lack of near infrared flux must be due to a large opacity towards the central star, contrary to other post-AGB objects where the lack of near infrared flux is due to the lack of dust close to the star (e.g. Hoogzaad et al. 2002). The high mass-loss rate suggested by our preliminary model is consistent with this large opacity, and can thus explain the lack of near infrared flux. However, more detailed modelling in the form of 2D radiative transfer calculations is needed in order to properly constrain the mass-loss rate of IRAS 16342.

The high mass-loss discussed above is only expected to occur during the SW phase of the AGB. Combined with the proto-PN nature of IRAS 16342, we conclude that IRAS 16342 must be a very young proto-PN, perhaps the youngest one observed so far. The star seems to be right at the transition point between the AGB phase (which explains its extreme OH/IR type spectrum) and the proto-PN phase. This makes IRAS 16342 a very interesting object for the study of the transition between OH/IR type AGB stars and proto-PN or post-AGB stars. The young proto-PN nature of IRAS 16342 is in agreement with SAH99. Still, the possibility of a more evolved central star can not be excluded, since the $`\mathrm{H}_2\mathrm{O}`$ maser velocities in the lobes (see Sect. 1) are more in favour of an A or late type B star. It is difficult to understand how a later type star can drive such a fast wind.

## 4 The infrared images

### 4.1 Description of the images

In Fig. 3 we present the deconvolved ISAAC and TIMMI2 images of IRAS 16342. The images were taken in the L, N7.9, N9.8, N2 and Q bands. The central wavelengths of these bands, $`\lambda _0`$, are given in Table 3. The lower right image shows the division of the N2 band image by the Q band image.

IRAS 16342 reveals a bipolar morphology in the L band. The bipolar nature is similar to that seen in the optical images taken with the HST by SAH99. In both cases we see two lobes, surrounded by a faint halo that follows the contours of the lobes. The TIMMI2 images reveal an ellipsoidal morphology, which may represent the bipolar nebula in the ISAAC and HST images, but at lower spatial resolution. Some caution is needed here however. The optical seeing during the observations with ISAAC and TIMMI2 were $`1.22.0^{\prime \prime }`$ and $`0.92.0^{\prime \prime }`$ respectively (see Sect. 2.2 and Sect. 2.3). In general, the seeing is better at longer wavelengths under the same optical seeing. Therefore we might also expect the TIMMI2 images to have an equal or even better spatial resolution than the ISAAC image, in which case differences seen between the images are intrinsic to IRAS 16342.

For each image, the size $`d`$ of the nebula along its major axis is listed in Table 3, together with $`P.A.`$, its position angle. We define the position angle as the angle between the major axis of the bipolar nebula and the north-south line, where $`0^{}`$ represents north and $`90^{}`$ represents east. Contrary to the size of the nebula, $`P.A.`$ changes significantly as a function of wavelength, showing a decrease with increasing wavelength. For the TIMMI2 images, both $`d`$ and $`P.A.`$ were measured by fitting two dimensional Gaussians to the images. The Full Width at Half Maximum (FWHM) of the Gaussian along its major axis and its tilt angle (i.e. $`P.A.90^{}`$), were used to obtain d and $`P.A.`$ respectively. For the ISAAC image we estimated $`d`$ and $`P.A.`$ by eye, while for the HST image we either used values given by SAH99 or derived these values from their data.

The N7.9 band and the Q band reveal (sub)structure additional to the ellipsoidal structure mentioned in the previous paragraph. The N7.9 band shows some extra extended emission in the north-east direction. Also, compared to the other bands, the morphology in the N7.9 band seems to be less smooth, with more spatial brightness variations. However, due to the low signal-to-noise quality of the data it is difficult to determine if the structures are real. Therefore, they should be treated with caution. The Q band image reveals extended emission, or a halo, surrounding the central ellips. The halo is roughly spherically symmetric, with a gap in the north eastern direction, and approximately 6′′ in diameter, i.e. twice the size of the bipolar nebula. The halo can also be seen in the lower right part of Fig. 3, displaying the division of the N2 band and the Q band images. The division shows that the halo is not uniformly bright.

In general our infrared images nicely agree with those made by Meixner et al. (1999) (also see Sect. 1). The halo we claim to see at $`20.0\mu `$m has not been reported by Meixner et al. (1999) however. It will therefore be interesting to make additional observations at $`20.0\mu `$m in the future.

### 4.2 The morphology of IRAS 16342

#### 4.2.1 The lobes

Given the similarities with the HST image, it seems natural to explain the bipolar nebula in the L band as due to reflection as well. Using this assumption, we will show in Sect. 4.3 that the dust particles in the lobes must be small ($`0.1\mu `$m). If scattering by small particles is indeed responsible for the presence of the bipolar nebula in the ISAAC image at $`3.78\mu `$m, a very red bright source is needed to illuminate the nebula, since small dust particles usually scatter poorly at wavelengths large compared to their radii. This very red source may either be the central star of IRAS 16342 or the inner regions of the torus. If it is the central star, this suggests that IRAS 16342 only recently entered the post-AGB phase, in agreement with our results from Sect. 3.3.

An alternative to the scattering scenario in the L band might be that the observed light is caused by thermal dust emission. This is unlikely however. The dust temperature in the lobes at a distance $`r`$ from the central star, $`T_\mathrm{d}`$, may be represented by $`T_\mathrm{d}=T_{\mathrm{eff}}(r/R_{})^{0.4}`$. This expression assumes that the lobes are optically thin for photospheric radiation and that the dust emission efficiency scales as $`Q(\lambda )\lambda ^1`$. Using Wiens law, we need $`T_\mathrm{d}\mathrm{\hspace{0.17em}650}\mathrm{K}`$ to emit efficiently at $`3.78\mu `$m. With $`T_{\mathrm{eff}}\mathrm{\hspace{0.17em}2670}\mathrm{K}`$ and $`R_{}=372\mathrm{R}_{}`$ (also see Sect. 3.3) we find that the radius of the region that can be heated to this temperature is $`r\mathrm{\hspace{0.17em}60}\mathrm{AU}`$. Using the distance of 2 kpc given by SAH99 and an angular diameter of 3′′, the observed radius of the lobes is 3000 AU. This is 50 times larger than the expected typical size of the emitting region derived above. We therefore exclude thermal dust emission as the source of radiation in the ISAAC image. Non equilibrium emission from small grains is also unlikely, since this requires the presence of UV photons (Andriesse 1978) and these are not likely to be present in the system. However, it can not be ruled out completely.

The behaviour of $`P.A.`$ as a function of wavelength suggests that (as expected) the dominant source of radiation changes from scattering to thermal emission. At the shortest wavelengths the source is closely aligned in the direction of the bipolar scattering nebulae seen in the HST images; at longer wavelengths the source tends to align more with the torus described in Sect. 1, which is expected to be dominated by thermal emission. Still, even up to $`20\mu `$m the elongation is approximately along the major axis of the nebula seen with HST and ISAAC, suggesting that the torus must be very cold ($`<150\mathrm{K}`$). This is in agreement with the very red ISO spectrum, which is most likely dominated by thermal emission from the material in the torus.

#### 4.2.2 The halo

We expect the halo in the Q band image to be due to a previous, more spherical, episode of mass-loss on the AGB. Given the size of the halo (6′′) and the bipolar nebula (3′′) and assuming an outflow velocity of $`15\mathrm{km}\mathrm{s}^1`$ during the AGB phase, and a distance of 2 kpc (SAH99), we find that the mass in the halo was ejected between approximately 950 and 1900 years ago. The Q band image therefore suggests that the transition from a (more) spherical to a more axial symmetric form of mass-loss for IRAS 16342 only occured near the very end of the AGB. We speculate that the non uniform brightness distribution may indicate that there are density enhancements in the halo, but this remains questionable. More observations are required to investigate this possibility.

### 4.3 Limits on grain sizes

Using a simple model, we can put constraints on the particle size of the dust in both the torus and the lobes. Now we will explain the model in detail and discuss its results. We will show that there are two independent methods that allow us to derive the optical depth in the torus along the line-of-sight towards the eastern lobe. The first method depends on the dust particle size in the lobes, while the second method depends on the particle size in the torus. Both methods must yield the same optical depth for a given wavelength. Applying this constraint and calculating the optical depth with both methods for three different wavelengths, yields the particle size in both the lobes and the torus.

Our model makes use of the following assumptions. The eastern lobe in the HST and ISAAC images is significantly redder than the western lobe. In order to explain this additional reddening, we assume the eastern lobe to be located completely behind the torus, and the western lobe to be positioned completely in front of the torus (see Fig. 4). The torus acts as a screen that causes additional reddening of the eastern lobe with respect to the western lobe. We assume that scattering is the dominant source of emission in the lobes, and that both scattering lobes are intrinsically equally bright. Even in the absence of the (extinction due to the) torus, the lobes will not appear equally bright as seen from Earth, due to the orientation of the system. The orientation causes the eastern lobe to be dominated by back scattering, while the western lobe is dominated by forward scattering. As a result, the eastern lobe will appear less bright than the western lobe. Interstellar reddening is assumed to be the same for both lobes, and scattering of stellar radiation within the torus is assumed to be unimportant (due to its large radial optical depth).

In the first method to derive the optical depth in the torus we measure the intensities of the eastern and western lobes, $`I_{e,\lambda }`$ and $`I_{w,\lambda }`$, as seen from Earth (see Fig. 4). These are given by

$$I_{e,\lambda }=I_0G_\lambda (180i,a_\mathrm{l})e^{\tau _\lambda }$$
(1)
$$I_{w,\lambda }=I_0G_\lambda (i,a_\mathrm{l})$$
(2)
respectively, where $`I_0`$ is the intrinsic intensity of each lobe. $`a_\mathrm{l}`$ is the characteristic particle size in the lobes, while $`G_\lambda (\theta ,a_\mathrm{l})`$ is the scattering phase function of a particle with radius $`a_\mathrm{l}`$ irradiated at wavelength $`\lambda `$ and scattering at an angle $`\theta `$. The phase function was calculated for size distributions of spherical grains using Mie theory (Bohren & Huffman 1983). Finally, $`i`$ is the inclination angle of the system, defined as the angle between the major axis of the nebula and the line-of-sight. We take $`i=40^{}`$ (SAH99). $`\tau _\lambda `$ is the optical depth of the torus along the line-of-sight towards the eastern lobe. From Eq. (1) and (2) we find

$$\tau _\lambda =\mathrm{ln}\left(\frac{I_{w,\lambda }G_\lambda (180i,a_\mathrm{l})}{I_{e,\lambda }G_\lambda (i,a_\mathrm{l})}\right)$$
(3)
This yields $`\tau _\lambda `$ for different values of $`a_\mathrm{l}`$.

The intensity ratio $`I_{w,\lambda }/I_{e,\lambda }`$ was determined from the magnitude difference between the lobes, $`\mathrm{\Delta }m`$, which was either given by SAH99 or we measured it ourselves. This implies that we replaced the intensity ratios by flux ratios. The values for the magnitude differences are listed in Table 4, together with the corresponding intensity ratios.

We describe the grain size distribution used in the calculation of $`G_\lambda (\theta ,a_\mathrm{l})`$ by a power law $`n(a)a^m`$ with $`m=3`$. This value was chosen for calculation conveniences. Differences with models using the more generally adopted ISM value of -3.5 (e.g. Mathis et al. 1977) are expected to be small. The effective width of the size distribution was kept fixed. The grains in the lobes were assumed to be composed of silicates, for which we used the optical constants of ‘cold’ silicates derived by Suh (1999). Using the optical constants of ‘warm’ silicates derived by Suh (1999) yields very similar results (not shown).

In the second method we derive $`\tau _\lambda `$ from the well known relation

$$\tau _\lambda =n_d\pi a_\mathrm{t}^2Q_\lambda ^{\mathrm{ext}}(a_\mathrm{t})L$$
(4)
where $`n_d`$ is the dust particle number density in a dust column of length $`L`$ along the line-of-sight; $`a_\mathrm{t}`$ and $`Q_\lambda ^{\mathrm{ext}}(a_\mathrm{t})`$ are the characteristic radius and extinction efficiency of the particles in the torus. For the calculation of the extinction efficiency we adopt the same dust properties as those assumed above for the calculation of the scattering phase functions.

Eq. (3) and (4) should yield the same result for $`\tau _\lambda `$ at a given wavelength. Ratios between optical depths are independent of $`n_d`$ and $`L`$ and therefore allow a relatively easy derivation of the particle sizes in the lobes and the torus.

Using Eq. (3) and the appropriate intensity ratios (see Table 4) and phase functions, we measured $`\tau _\lambda `$ at $`0.55`$, $`0.79`$ and $`3.78\mu `$m for different assumptions for the characteristic particle size in the lobes, $`a_\mathrm{l}`$. In Fig. 5 $`\frac{\tau _{0.55}}{\tau _{0.80}}`$ is plotted versus $`\frac{\tau _{0.80}}{\tau _{3.78}}`$ with $`a_\mathrm{l}`$ as a free parameter (see line indicated with ‘I-ratios’). We also calculated $`Q_\lambda ^{\mathrm{ext}}(a_\mathrm{t})`$ at the above wavelengths for different values of the characteristic particle size in the torus, $`a_\mathrm{t}`$. Since $`\tau _\lambda Q_\lambda ^{\mathrm{ext}}(a_\mathrm{t})`$, we must have

$$\frac{\tau _{0.55}}{\tau _{0.80}}=\frac{Q_{0.55}^{\mathrm{ext}}(a_\mathrm{t})}{Q_{0.80}^{\mathrm{ext}}(a_\mathrm{t})}$$
(5)
and

$$\frac{\tau _{0.80}}{\tau _{3.78}}=\frac{Q_{0.80}^{\mathrm{ext}}(a_\mathrm{t})}{Q_{3.78}^{\mathrm{ext}}(a_\mathrm{t})}$$
(6)
independent of $`n_d`$ and $`L`$. In Fig. 5 we again plot $`\frac{\tau _{0.55}}{\tau _{0.80}}`$ versus $`\frac{\tau _{0.80}}{\tau _{3.78}}`$, this time with $`a_\mathrm{t}`$ as a free parameter (see line indicated with ‘Q-ratios’). The intersection of the two curves provides a consistent set for $`a_\mathrm{l}`$ and $`a_\mathrm{t}`$.

For the lobes we find that we have a size distribution with a maximum particle size $`a_\mathrm{l}^{\mathrm{max}}=0.09\pm 0.02\mu `$m, while for the torus we find $`a_\mathrm{t}^{\mathrm{max}}=1.3\pm 0.2\mu `$m. Both values are independent of the adopted effective width of the size distribution, which we take to be the same for particles in both the torus and the lobes. A test with amorphous olivines instead of the cold ‘Suh’ silicates yielded a similar result. The errors on the particle sizes are derived from an error analysis on Eq. (3) taking the errors on $`\mathrm{\Delta }m`$ and $`I_{w,\lambda }/I_{e,\lambda }`$ into account for each wavelength (see Table 4). We are interested in the maximum particle size of the distribution instead of the typical (i.e. effective) particle size. This is because we want to set upper limits on the particle sizes found in the lobes and torus. The typical size simply scales linearly with the maximum particle size.

We thus derive that the maximum particle size in the torus is about an order of magnitude larger than the maximum particle size in the lobes. In an AGB outflow, the particle size of the grains formed in the outflow is proportional to the density of the envelope at the dust condensation radius (Dominik 1992). Larger densities provide larger grains. The spatial distribution of large particles in a torus geometry suggests that this material is expelled from the AGB star in an aspherical high density wind. In a natural way, this accounts for the presence of large grains in the torus and small grains in the lobes. However, the bipolar morphology of IRAS 16342 may also have originated from a later phase, i.e. after the AGB. In this case, the AGB wind may have been spherical, but due to the interaction of this wind with the (more recent) fast bipolar molecular outflow, the morphology may have changed into the observed torus like spatial distribution. The interaction of these two flows may have caused a shattering of the dust particles in the lobes into smaller ones, while leaving the grains in the torus intact. Regardless of the correct scenario, the presence of the bipolar nebula in IRAS 16342 at optical wavelengths suggests that the mass-loss rate near the central star very recently must have dropped, otherwise the optical photons in this nebula would not be able to escape. Obviously the star has already managed to create a bipolar cavity while the optical depth of the dust shell in the equatorial torus is still huge.

The derived particle size for the torus is large compared to the typical ISM values, which are of the order of $`a\mathrm{\hspace{0.17em}0.01}`$ to $`\mathrm{\hspace{0.17em}0.2}\mu `$m (e.g. Mathis et al. 1977), but in agreement with other objects for which large mass-loss rates are expected. Examples of such sources are the binary system AFGL 4106 (Molster et al. 1999) and the post-AGB star HD161796 (Hoogzaad et al. 2002). The large grain size in the torus explains why the dust particles in the torus of IRAS 16342 are grey compared to the ISM.

We finally note that when the emission in the ISAAC image at $`3.78\mu `$m is caused by non equilibrium emission from small dust grains, as was considered in Sect. 4.2.1, the resulting particle sizes in the lobes and torus are very similar to those derived above. The possibility of non equilibrium emission by small grains at $`3.78\mu `$m can simply be modeled by leaving $`F_{3.78}(180i,a_\mathrm{l})`$ and $`F_{3.78}(i,a_\mathrm{l})`$ out of equation 3 for $`\tau _{3.78}`$, i.e. by assuming that at $`3.78\mu `$m there is no scattering.

## 5 Conclusions

The results of our study can be summarised as follows:

1. The ISO spectrum of IRAS 16342 is that of an extremely reddened OH/IR star. In addition to amorphous silicate absorption features at $`10`$ and $`20\mu `$m, the spectrum shows crystalline silicate features in absorption at wavelengths up to $`45\mu `$m. Only at longer wavelengths the crystalline features are in emission. No other OH/IR star is known to show this behaviour. We estimate the mass-loss rate of IRAS 16342 may be as large as $`10^3\mathrm{M}_{}\mathrm{yr}^1`$. This high rate explains the lack of flux in the near-IR as well as the peak wavelength in the SED.
2. We identify forsterite ($`\mathrm{Mg}_2\mathrm{SiO}_4`$), diopside ($`\mathrm{MgCaSi}_2\mathrm{O}_6`$), and possibly clino-enstatite ($`\mathrm{MgSiO}_3`$) in the spectrum. Besides these crystalline silicates, also crystalline water ice is detected. All these materials are also found around other (high $`\dot{M}`$), oxygen rich, evolved stars. At $`48\mu `$m a (to date) unidentified feature is seen.
3. IRAS 16342 must be a very young proto-PN, perhaps the youngest one observed until present. The central star seems to be right at the transition point between the AGB phase and the proto-PN phase.
4. The deconvolved L band ISAAC image of IRAS 16342 at $`3.78\mu `$m reveals a bipolar scattering nebula. Both lobes of the nebula are separated by a dark equatorial waist, or torus. This morphology is consistent with that found in HST images. The presence of the scattering nebula in the L band suggests that the central star of IRAS 16342 may be very red, which would be in agreement with the young proto-PN nature of the system. The TIMMI2 Q band image reveals thermal emission from an extended, roughly spherical halo of $`\mathrm{\hspace{0.17em}6}^{\prime \prime }`$ in diameter. We expect the halo to be due to a previous, more spherical, phase of mass-loss on the AGB. This halo material was ejected between approximately 950 and 1900 years ago.
5. TIMMI2 images at wavelengths intermediate to the L and Q band reveal an ellipsoidal morphology, which may represent the bipolar nebula but at lower spatial resolution. The position angle of the nebula changes as a function of wavelength, suggesting that the nature of the radiation changes from scattered light to thermal emission (by cold material). Still, even up to $`20\mu `$m the elongation is approximately along the major axis of the nebula seen in the optical and near-IR, suggesting that the equatorial torus must be very cold. This is in agreement with the very red ISO spectrum, which is most likely dominated by thermal emission from this torus.
6. Using a simple model for the geometry of IRAS 16342, we calculate the dust particle sizes in both the torus and the lobes. The grains in the lobes have a maximum particle size $`a_\mathrm{l}^{\mathrm{max}}=0.09\mu `$m, while for the torus we find $`a_\mathrm{t}^{\mathrm{max}}=1.3\mu `$m. The relatively large torus grains suggest that the bipolar morphology of IRAS 16342 may already have been created during the AGB phase.

###### Acknowledgements.

C. Dijkstra wishes to thank Sacha Hony for his help on the data reduction of the TIMMI2 data. We are grateful for the support by the staff of the ESO 3.6m telescope.


format: text

text:
# Modified two-potential approach to tunneling problems

## I Introduction

The quantum mechanical tunneling through a classically forbidden region is an ubiquitous phenomenon in physics, which has been extensively studied since the early days of quantum mechanics. In 1927, Hund Hund  was the first to point out the possibility of “barrier penetration” between two discrete states. In the same year, Nordheim Nordheim  considered the case of tunneling between continuum states. Subsequently, Oppenheimer Oppy  performed a calculation of the rate of ionization of the hydrogen atom, and Gamow, Gamow  Gurney and Condon Gurney  explained alpha decay rates of radioactive nuclei in terms of the tunneling effect.

While the semi-classical treatment of tunneling turned out to be very successful in many applications, the numerical calculation offers very little insight into the physical process. In addition, the validity of the standard WKB formula is rather restricted. Other methods, although more accurate, contain various uncertainties. For example, the results of the commonly used R-matrix theory lane  are often sensitive to the choice of the matching radius humb ; Arima , and the theoretical error is difficult to estimate.

The treatment of the tunneling problem can be essentially simplified by reducing it to two separate problems: a bound state problem and a non-resonant (scattering) state problem. This can be done consistently in the two-potential approach (TPA) gk ; g ; g1  (see also Refs. jrb ; tore ), representing the barrier potential as a sum of the “inner” and the “outer” terms, containing only bound and only scattering states, respectively. This approach not only provides better physical insights than many other approximations but it is also simple and accurate.

In this paper we propose further developments and a modification of the TPA, and present a detailed comparison of this approach with the results of numerical calculations based on the Gamow-state (resonant-state) formalism. The resulting analytical expressions are easy to interpret and they can be straightforwardly extended to the non-spherical case.

The paper is organized as follows. In Sect. II, the TPA is briefly described. Section III deals with the quantal correction terms to the TPA. The minimization of these terms prescribes unambiguously the “window” for the separation radius that divides the original barrier potential into inner and outer terms. In this case, by considering examples of wide and narrow nuclear resonances, we demonstrate that the TPA yields results which are practically the same as those of the resonant-state calculation. In Sec. IV, we present a modification of the TPA. The resulting expressions resemble those of the R-matrix theory, yet without any uncertainties related to the matching radius, see Sec. V. Finally, the summary of our work is contained in Sec. VI.

## II Two-potential approach

Consider a quantum well $`V(r)`$ with a barrier, which contains a quasi-stationary state at the $`E_{res}`$.

The coordinate space can be divided into two regions, the “inner” region, $`0<r<R`$, and the “outer” region, $`r`$$`>`$$`R`$, where $`R`$ is taken inside the barrier (see Fig. 1). Accordingly, one can introduce the two auxiliary potentials: the inner potential

$$U(r)=\{\begin{array}{cc}V(r)\hfill & \text{for}rR\hfill \\ V(R)=V_0\hfill & \text{for}r>R\hfill \end{array}$$
(1)
and the outer potential

$$\stackrel{~}{W}(r)=\{\begin{array}{cc}V_0\hfill & \text{for}rR\hfill \\ V(r)\hfill & \text{for}r>R\hfill \end{array}.$$
(2)
The inner potential contains a bound state, $`\mathrm{\Phi }_0(r)`$ ($`E`$=$`E_0`$), representing an eigenstate of the “inner” Hamiltonian $`H_0=K+U(r)`$, where $`K=^2/2m`$ is the kinetic energy term ($`\mathrm{}`$=1). One can demonstrate g  that the energy and the width of the quasi-stationary state, associated with the complex-energy poles of the total Green’s function $`G(E)=(EKV)^1`$, are obtained from the following equation

$$E=E_0+\mathrm{\Phi }_0|W|\mathrm{\Phi }_0+\mathrm{\Phi }_0|W\stackrel{~}{G}(E)W|\mathrm{\Phi }_0.$$
(3)
Here $`W(r)=\stackrel{~}{W}(r)V_0`$, and the Green’s function $`\stackrel{~}{G}`$ is given by

$$\stackrel{~}{G}(E)=G_0(E)\left[1+\stackrel{~}{W}\stackrel{~}{G}(E)\right],$$
(4)
where

$$G_0(E)=\frac{1\mathrm{\Lambda }}{E+U_0KU},\mathrm{\Lambda }=|\mathrm{\Phi }_0\mathrm{\Phi }_0|.$$
(5)
The resonance energy $`E_{\mathrm{res}}=\mathrm{}(E)`$ and the width $`\mathrm{\Gamma }=2\mathrm{}(E)`$ of a quasi-stationary state obtained from Eq. (3) are independent of the choice of the separation radius $`R`$.

Equation (3) can be solved iteratively by using the standard Born series for the Green’s function $`\stackrel{~}{G}`$, i.e. by expanding $`\stackrel{~}{G}`$ in powers of $`G_0`$. Yet, the corresponding expansion for the quasi-stationary state energy converges very slowly. For that reason, we proposed g  a more efficient expansion scheme in which $`\stackrel{~}{G}`$ is expanded in powers of the Green’s function $`G_{\stackrel{~}{W}}(E)=(EK\stackrel{~}{W})^1`$, corresponding to the outer potential $`\stackrel{~}{W}`$. From Eq. (4) it immediately follows that

$$\stackrel{~}{G}=G_{\stackrel{~}{W}}+G_{\stackrel{~}{W}}\left(UU_0\right)\stackrel{~}{G}\stackrel{~}{G}_{\stackrel{~}{W}}\mathrm{\Lambda }\left(1+\stackrel{~}{W}\stackrel{~}{G}\right),$$
(6)
Iterating Eq. (6) in powers of $`G_{\stackrel{~}{W}}`$ and then substituting the result into (3), one finds the desirable perturbative expansion for the energy and the width of the resonance. By truncating this series, one obtains the following first-order relation valid for the isolated metastable state:

$$E=E_0+\mathrm{\Phi }_0|W|\mathrm{\Phi }_0+\mathrm{\Phi }_0|WG_{\stackrel{~}{W}}(E)W|\mathrm{\Phi }_0.$$
(7)
The above equation can be solved iteratively for $`E=E_ri\mathrm{\Gamma }/2`$ by assuming that the energy shift $`\mathrm{\Delta }=E_rE_0`$ and the width $`\mathrm{\Gamma }`$ are small compared to $`E_0`$ and $`V_0E_0`$. In such a case, one can put $`G_{\stackrel{~}{W}}(E)G_{\stackrel{~}{W}}(E_0)`$, thus reducing Eq. (7) to

$$E=E_0+\mathrm{\Phi }_0|W|\mathrm{\Phi }_0+\mathrm{\Phi }_0|WG_{\stackrel{~}{W}}(E_0)W|\mathrm{\Phi }_0.$$
(8)
By using the Schrödinger equation for $`G_{\stackrel{~}{W}}`$ one finally obtains the TPA expressions gk ; g  for the width $`\mathrm{\Gamma }`$ and for the energy shift $`\mathrm{\Delta }=E_{\mathrm{res}}E_0`$ of the quasi-stationary state:

$`\mathrm{\Gamma }`$ $`=`$ $`{\displaystyle \frac{1}{mk}}\left[\mathrm{\Phi }_0(R)\chi _k^{}(R)\mathrm{\Phi }_0^{}(R)\chi _k(R)\right]^2,`$ (9)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }_0^2(R)}{2mk}}\left[\alpha \chi _k(R)+\chi _k^{}(R)\right]\left[\alpha \stackrel{~}{\chi }_k(R)+\stackrel{~}{\chi }_k^{}(R)\right],`$ (10)
where $`k=\sqrt{2mE_0}`$, $`\alpha =\sqrt{2m(V(R)E_0)}`$, $`\stackrel{~}{\chi }_k(r)=\mathrm{}(\chi _k^{(+)}(r))`$, and $`\chi _k`$ ($`\chi _k^{(+)}`$) stands for the irregular (outgoing) solution of the Schrödinger equation for the outer potential.

It follows from Eqs. (9) and (10) that both $`\mathrm{\Gamma }`$ and $`\mathrm{\Delta }`$ are given in terms of bound and scattering state wave functions. Thus TPA essentially simplifies the treatment of tunneling, because the standard approximation schemes can be used for evaluation of $`\mathrm{\Phi }_0`$ and $`\chi _k`$. For instance, by applying the semi-classical approximation, one obtains the improved Gamow formula for $`\mathrm{\Gamma }`$ gk ; g , which is useful for different applications buck ; nazar . In particular, an extension of Eqs. (9) and (10) to the multi-dimensional case can be found in Ref. g1 .

## III Corrections to TPA and the choice of the separation radius

The accuracy of Eqs. (9) and (10) can be determined by evaluating the leading correction terms. There are two types of corrections to TPA: (a) those due to the replacement of $`\stackrel{~}{G}`$ by $`G_{\stackrel{~}{W}}`$ in Eq. (3) leading to Eq. (7), and (b) those due to the replacement of $`G_{\stackrel{~}{W}}(E)`$ by $`G_{\stackrel{~}{W}}(E_0)`$ in Eq. (7) leading to Eqs. (9), (10). The correction terms of the first type, $`(\mathrm{\Delta }\mathrm{\Gamma })_1`$, can be obtained by iterating (6). One finds from the first iteration g :

$$\frac{(\mathrm{\Delta }\mathrm{\Gamma })_1}{\mathrm{\Gamma }}\frac{mV^{}(R)}{16[2m(V_0E_0)]^{3/2}}.$$
(11)
Equation (11) might suggest that an optimal choice of the separation radius corresponds to $`V^{}(R)=0`$, i.e., the maximum of $`V(r)`$. However, it has been demonstrated numerically nazar ; talou  that if the top of the barrier is close to the closing potential, such a choice is not optimal, since in this case Eqs. (9) and (10) become less accurate. The reason is that the energy shift $`\mathrm{\Delta }`$ becomes appreciable so that (7) cannot be replaced by (8).

This can be illustrated by considering a square-well potential discussed in Ref. erich : $`V(r)=l(l+1)/2mr^2U_0`$ for $`r<R_1`$, and $`V(r)=l(l+1)/2mr^2`$ for $`rR_1`$, where the top of the barrier coincides with the closing potential (see Fig. 2).

For a P-wave resonance, the numerical calculation gives $`E_{\mathrm{res}}=1`$ keV and $`\mathrm{\Gamma }=55.9`$ eV erich . Now we apply the TPA by taking the separation radius at the boundary, $`R=R_1`$. The corresponding inner potential $`U(r)`$ has a bound state at $`E_0`$=275 keV. However, the corresponding energy shift, $`\mathrm{\Delta }`$=$``$300 keV, is of the same order of magnitude as the energy $`E_0`$ of the bound state. Consequently, the replacement of $`G_{\stackrel{~}{W}}(E)`$ by $`G_{\stackrel{~}{W}}(E_0)`$ in (7) leads to large corrections to the resonance energy and the width so that Eqs. (9) and (10) cannot be used.

Let us estimate the correction term $`(\mathrm{\Delta }\mathrm{\Gamma })_2=\mathrm{\Gamma }(E_0+\mathrm{\Delta })\mathrm{\Gamma }(E_0)`$ due to such a replacement. One can use the semi-classical Gamow formula, $`\mathrm{\Gamma }\mathrm{exp}(2_{r_1}^{r_2}|p(r)|dr`$, with $`r_1`$ and $`r_2`$ being the inner and outer classical turning points, respectively, and $`|p(r)|=\sqrt{2m[V(r)E_0]}`$. Approximating $`V(r)`$ for $`r_1<r<r_2`$ by the inverted harmonic oscillator, one obtains

$$\frac{(\mathrm{\Delta }\mathrm{\Gamma })_2}{\mathrm{\Gamma }}\mathrm{\Delta }_{r_1}^{r_2}\frac{\sqrt{2m}}{\sqrt{V(r)E_0}}𝑑r\pi \frac{\sqrt{2m}(r_2r_1)}{2\sqrt{\overline{V}E_0}}\mathrm{\Delta },$$
(12)
where $`\overline{V}`$= max $`V(r)`$ and $`\mathrm{\Delta }`$ is is given by Eq. (10).

Thus, in order to reduce the correction term $`(\mathrm{\Delta }\mathrm{\Gamma })_2`$, one needs to minimize the energy shift $`\mathrm{\Delta }`$. It follows from (10) that $`\mathrm{\Delta }`$ can be strongly suppressed by taking the separation radius $`R`$ deeply inside the barrier. Indeed, Eq. (10) contains a product of regular and irregular wave functions, which do not vary considerably under the barrier. However, the factor $`|\mathrm{\Phi }_0(R)|^2`$ decays exponentially with $`R`$. Therefore, by taking the separation radius $`R`$ far away from the boundary, $`Rr_1`$, one finds that $`\mathrm{\Delta }0`$, and $`E_0E_{\mathrm{res}}`$. As a result, (7) can be replaced by (8), leading to Eqs. (9) and (10).

To illustrate this point, let us again consider the example of a P-wave resonance discussed above. By taking the separation radius $`R>R_1`$, we readily find that $`\mathrm{\Delta }0`$ as $`RR_1`$ increases, and $`E_0E_{\mathrm{res}}`$= 1 keV. The resulting value of $`\mathrm{\Gamma }`$=55.3 eV is very close to the Gamow-state value of the resonance width. The separation radius $`R`$ cannot be chosen too close to the outer classical turning point since in such a case $`V(R)E_0`$ and the correction term (11) becomes important. In fact, Eqs. (11) and (12) define the lower and upper limits of $`R`$.

In the following, we discuss the TPA results for proton and neutron resonances in the realistic average nuclear potential. The approximate TPA expressions are compared to the resonant states results obtained using the GAMOW code gamow .

### III.1 Comparison with Gamow-state calculations

Consider single-nucleon resonances in the Woods-Saxon (WS) potential $`V_{WS}(r)`$ represented by a sum of central, spin-orbit, centrifugal, and Coulomb terms. Here, we apply the parametrization of Ref. nazar , namely: $`R_0`$=1.17A1/3 fm, $`a=0.75`$ fm for the central term, and $`U_0^{so}=0.2U_0`$ and $`R_0^{so}=1.01`$A1/3 for the spin-orbit potential. We calculate $`\mathrm{\Delta }`$ and $`\mathrm{\Gamma }`$ according to Eqs. (9) and (10) by varying the separation radius $`R`$ inside the barrier, starting with the barrier radius, $`\overline{R}`$, corresponding to the maximum of $`V(r)`$.

We begin with the high-$`\mathrm{}`$ narrow proton resonance 0h11/2. The parameters of the WS potential are appropriate for 147Tm, which is a proton emitting nucleus. The potential depth, $`U_0`$=$``$61.8823 MeV, was adjusted to the energy $`E_{\mathrm{res}}`$=1.5 MeV. The resulting barrier radius is $`\overline{R}`$=8.54 fm ($`V(\overline{R})`$=17.44 MeV) and the inner and outer turning points are $`r_1`$=6.33 fm and $`r_2`$=71.15 fm, respectively. Since $`E_{\mathrm{res}}V(\overline{R})`$, the calculated 0h11/2 resonant state has a very small width, $`\mathrm{\Gamma }`$=4.918 10-18 MeV. The results of TPA calculations are shown in Table I for different values of $`R\overline{R}`$, together with the corresponding correction terms to TPA: $`(\mathrm{\Delta }\mathrm{\Gamma })_1/\mathrm{\Gamma }_{\mathrm{TPA}}`$ (11) and $`(\mathrm{\Delta }\mathrm{\Gamma })_2/\mathrm{\Gamma }_{\mathrm{TPA}}`$ (12). Table I also displays the actual accuracy of the TPA, $`(\mathrm{\Delta }\mathrm{\Gamma })/\mathrm{\Gamma }`$, where $`(\mathrm{\Delta }\mathrm{\Gamma })=\mathrm{\Gamma }\mathrm{\Gamma }_{\mathrm{TPA}}`$. Since $`\overline{R}>r_1`$, the energy shift $`\mathrm{\Delta }`$ is small for $`R=\overline{R}`$. Therefore, the results of TPA are in a good agreement with the resonant-state calculations already for $`R=\overline{R}`$. Next we consider the low-$`\mathrm{}`$, broader 2s1/2 resonance at $`E_{\mathrm{res}}`$=1.5 MeV, which is considerably closer to the top of the barrier $`V(\overline{R})=9.43`$ MeV ($`\overline{R}`$=9.34 fm). As shown in Table I, also in this case $`\mathrm{\Gamma }_{\mathrm{TPA}}`$ nicely agrees with the numerical result, and the accuracy of the TPA is well estimated by Eqs. (11) and (12).

Table II displays the TPA results for the 0i13/2 and 1f5/2 neutron resonances in 133Sn at an energy $`E_{\mathrm{res}}`$=1 MeV. Here $`\mathrm{\Gamma }`$ is much larger due to the absence of the Coulomb barrier. As in the proton case, there is very good agreement with numerical calculations, provided that $`R`$ is taken far away from the turning points, inside the window determined by (11) and (12), and the results of TPA weakly depend on the separation radius $`R`$. This suggests that the separation radius can be eliminated altogether from the TPA expressions. As demonstrated in the following section one can indeed modify the TPA in such a way.

## IV Modified Two-Potential Approach

A tunneling potential can always be written as a sum of attractive and repulsive parts, $`V(r)=V_{att}(r)+V_{rep}(r)`$, where $`V_{rep}(r)`$ becomes dominant at distances beyond the barrier radius (see Fig. 3).

Therefore, starting with some radius $`\overline{r}`$, the total potential $`V(r)`$ can be well approximated by the repulsive component only. The value of $`\overline{r}`$ should be chosen in such a way that the attractive (e.g., nuclear) part can disregarded with a desired accuracy $`\eta 1`$:

$$|1V_{rep}(r)/V(r)|\eta \text{for}r\overline{r}.$$
(13)
For instance, in the case of a square-well potential of Fig. 2, Eq. (13) is satisfied for any $`\eta `$ provided that $`\overline{r}`$$`>`$$`R_1`$. In most cases, the attractive potential rapidly decreases beyond the barrier radius, so that $`\overline{r}`$ is closer to $`\overline{R}`$ than to the separation radius $`R`$ in the TPA (cf. Fig. 3).

Usually, the repulsive part $`V_{rep}(r)`$ is well known, as well as the two linearly independent (regular and irregular) solutions $`F_k(r)`$ and $`G_k(r)`$ of the corresponding Schrödinger equation. For instance, if $`V_{rep}(r)`$ is a sum of Coulomb and centrifugal potentials, then $`F_k(r)`$ and $`G_k(r)`$ are the standard Coulomb functions. This implies that any solution of the Schrödinger equation with the potential $`V(r)`$ can be written for $`r>\overline{r}`$ as the linear combination of $`F_k(r)`$ and $`G_k(r)`$.

Consider the bound-state wave function $`\mathrm{\Phi }_0(r)`$ of the inner potential $`U(r)`$ of Eq. (1). Since $`U(r)V_{rep}(r)`$ for $`\overline{r}rR`$, $`\mathrm{\Phi }_0(r)`$ can be expanded in this region as

$$\mathrm{\Phi }_0(r)=c_1G_k(r)+c_2F_k(r),$$
(14)
where $`k=(2mE_0)^{1/2}`$. The coefficients $`c_{1,2}`$ and the energy $`E_0`$ are obtained from matching of the logarithmic derivatives at $`r=\overline{r}`$ and $`r=R`$:

$`{\displaystyle \frac{c_1G_k^{}(R)+c_2F_k^{}(R)}{c_1G_k(R)+c_2F_k(R)}}`$ $`=`$ $`|p(R)|,`$ (15a)
$`{\displaystyle \frac{c_1G_k^{}(\overline{r})+c_2F_k^{}(\overline{r})}{c_1G_k(\overline{r})+c_2F_k(\overline{r})}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }_0^{}(\overline{r})}{\mathrm{\Phi }_0(\overline{r})}}.`$ (15b)
Note that $`\mathrm{\Phi }_0(r)\mathrm{exp}[|p(R)|(rR)]`$, for $`r>R`$. Solving Eqs. (15) one easily finds

$$\frac{\mathrm{\Phi }_0^{}(\overline{r})}{\mathrm{\Phi }_0(\overline{r})}=\frac{G_k^{}(\overline{r})}{G_k(\overline{r})}\left[\frac{1+\alpha a_k^1f_1(1+\stackrel{~}{a}_k^1)}{1+\alpha a_k^1f_2(1+\stackrel{~}{a}_k^1)}\right],$$
(16)
where $`\alpha =|p(R)|`$, $`a_k=F_k^{}(R)/F_k(R)`$, $`\stackrel{~}{a}_k=G_k^{}(R)/G_k(R)`$, and

$$f_1=\frac{G_k^{}(R)F_k^{}(\overline{r})}{F_k^{}(R)G_k^{}(\overline{r})},f_2=\frac{G_k^{}(R)F_k(\overline{r})}{F_k^{}(R)G_k(\overline{r})}.$$
(17)
The wave functions $`F_k(r)`$ and $`G_k(r)`$ are of the same order of magnitude in the asymptotic region, $`rr_2`$. However, in the classically forbidden region the regular wave function exponentially decreases and the irregular one exponentially increases with decreasing $`r`$. Using the semi-classical approximation, one can estimate $`G_k(R)F_k(\overline{r})\mathrm{exp}[_{\overline{r}}^R|p(r)|𝑑r]`$. Therefore, the coefficients $`f_{1,2}`$ are of the order of $`\mathrm{exp}[\alpha (R\overline{r})]`$, so the corresponding terms $`f_{1,2}(1+\stackrel{~}{a}_k^{(1)})`$ in Eq. (16) are exponentially suppressed and can be neglected. As a result, the matching condition (16) can be written as

$$\frac{\mathrm{\Phi }_0^{}(\overline{r})}{\mathrm{\Phi }_0(\overline{r})}=\frac{G_k^{}(\overline{r})}{G_k(\overline{r})}.$$
(18)
The above equation constitutes the MPTA condition for the resonance energy $`E_{\mathrm{res}}=E_0`$. In contrast to Eq. (16), the relation (18) does not exhibit any explicit $`R`$-dependence. Consequently, there is no need to evaluate the bound-state wave function at large values of $`R`$ well inside the barrier. Note also that for a narrow resonance the irregular wave function $`G_k(r)`$ is proportional to the real part of the outgoing (Gamow) solution:

$$\mathrm{\Psi }_{k_{res}}^{out}(r)G_{k_{res}}(r)+iF_{k_{res}}(r)=O_{k_{res}}(r),$$
(19)
with complex $`k_{res}=ki\gamma `$ and $`k_{res}^2=2m(E_{res}i\mathrm{\Gamma }/2)`$. Since for small $`\mathrm{\Gamma }`$ and $`\gamma `$ the imaginary parts of $`k_{res}`$ and $`G_{k_{res}}(r)`$ can be neglected, $`k_{res}k`$ and $`G_{k_{res}}(r)G_k(r)`$ (also for the outgoing Coulomb wave function $`O_{k_{res}}(r)=O_k(r)`$). In this case Eq. (18) represents the matching condition for the inner (bound state) wave function with the real part of the non-normalized outgoing wave. If the imaginary parts are not negligible, Eq. (18) can be still interpreted in terms of the standing-wave boundary condition at $`\overline{r}`$, which means that the scattering phase shift is $`\pi /2`$. The requirement that the phase shift is equal to $`\pi /2`$ represents an alternative definition for the position of a resonance in the absence of a non-resonant phase shift.

Consider now Eq. (9) for the width. Since $`R>\overline{r}`$, the outer wave function $`\chi _k(r)`$ in the region $`rR`$ can be represented by the linear combination of the regular and irregular solutions of $`V_{rep}(r)`$, and the corresponding coefficients in the linear combination of $`F_k(r)`$ and $`G_k(r)`$ are directly related to the (non-resonant) scattering phase shift for the outer potential $`\stackrel{~}{W}(r)`$. One easily finds

$$\chi _k(r)=\mathrm{cos}\delta _kF_k(r)+\mathrm{sin}\delta _kG_k(r)\text{for}rR,$$
(20)
where the phase shift $`\delta _k`$ is obtained from matching of logarithmic derivatives at the separation radius $`R`$:

$$\mathrm{tan}\delta _k=\frac{F_k(R)}{G_k(R)}\left(\frac{a_k\alpha }{\stackrel{~}{a}_k\alpha }\right).$$
(21)
Here we neglected the terms $`\mathrm{exp}(2\alpha R)`$. Substituting (14) and (20) into (9) and taking into account the Wronskian relation between $`F_k(r)`$ and $`G_k(r)`$ one obtains:

$$\mathrm{\Gamma }=\mathrm{cos}^2\delta _k\frac{k}{m}\left[\frac{\mathrm{\Phi }_0(\overline{r})}{G_k(\overline{r})}\right]^2\left[1\frac{(\stackrel{~}{a}_k+\alpha )(a_k\alpha )}{(a_k+\alpha )(\stackrel{~}{a}_k\alpha )}\right]^2.$$
(22)
Note that in the classically forbidden region $`F_k^{}(r)/F_k(r)|p(r)|`$ and $`G_k^{}(r)/G_k(r)|p(r)|`$, so that the second term in brackets of (22) can be neglected. In addition, $`\mathrm{cos}\delta _k1`$, as follows from (21). As a result, one arrives at the following simple expression for the width:

$$\mathrm{\Gamma }=\frac{k}{m}\left[\frac{\mathrm{\Phi }_0(\overline{r})}{G_k(\overline{r})}\right]^2.$$
(23)
Thus, similar to Eq. (18) for the resonance’s energy, the separation radius $`R`$ does not appear explicitly in the expression for the width.

Equations (18) and (23) represent the final result of the modified two-potential approach (MTPA). Despite their simple appearance, these expressions are very accurate. In fact, the accuracy of the MTPA is practically the same as that of the TPA since the former was derived from the latter by neglecting only small correction terms of the order of the accuracy of the TPA itself. For instance, for the previously discussed case of the P-wave resonance in the square well potential, one finds $`E_{\mathrm{res}}`$=1 keV and $`\mathrm{\Gamma }`$=55.3 eV, i.e. the same result as in TPA. In general, the accuracy of the MTPA can be estimated by means of the parameter $`\eta `$, which defines the lower limit for the matching radius $`\overline{r}`$. However, one has to keep in mind that $`\overline{r}`$ cannot be very large since the derivation of Eqs. (18) and (23) is valid only for $`\alpha (R\overline{r})1`$. Therefore, the value of $`\eta `$ is restricted by Eq. (11) in which $`R`$ is replaced by $`\overline{r}`$.

It is worth noting that an expression similar to Eq. (23) was used in Refs. maglione ; ferreira  for calculating partial widths for proton emission. The corresponding formula that applies to the single channel case can be written as

$$\mathrm{\Gamma }(r)=\frac{k}{m}\left|\frac{\mathrm{}\left(\mathrm{\Psi }_{k_{res}}^{out}(r)\right)}{O_k(r)}\right|^2,$$
(24)
where $`r`$ is large. The $`r`$-dependence of $`\mathrm{\Gamma }(r)`$ is weak but it can be reduced if one takes a more appropriate expression:

$$\mathrm{\Gamma }(r)=\frac{k}{m}\left|\frac{\mathrm{\Psi }_{k_{res}}^{out}(r)}{O_{k_{res}}(r)}\right|^2.$$
(25)
Another expression for the width can be derived from the continuity relation for the resonant states humb ; barmore :

$$\mathrm{\Gamma }(r)=i\frac{1}{2m}\frac{\mathrm{\Psi }_{k_{res}}^{out}(r)\mathrm{\Psi }_{k_{res}}^{out}(r)\mathrm{\Psi }_{k_{res}}^{out}(r)\mathrm{\Psi }_{k_{res}}^{out}(r)}{_0^r|\mathrm{\Psi }_{k_{res}}^{out}(r^{})|^2𝑑r^{}}.$$
(26)
This form is completely independent of $`r`$ in a wider range barmore ; vertseproc  and furnishes a value which is equal to that coming from the imaginary part of the energy. However, for very narrow resonances, it is difficult to calculate the imaginary part of the energy with sufficient precision. The expression (23) derived in the MTPA replaces the Gamow wave function with the (real) bound-state wave function $`\mathrm{\Phi }_0`$. Finally, let us emphasize that while Eq.(24) resembles the MTPA expression, it is based on different approximations and boundary conditions. On the other hand there are close connections between the MTPA and the R-matrix theory, also employing real-energy eigenstates, see Sec. V.

### IV.1 MTPA: numerical examples

We present below in Table III the results of the MTPA for the widths (23) of resonances discussed in Tables I and II in the context of TPA. (Since the MTPA resonance energies are very close to the exact result, they are not displayed.) One finds that the MTPA reproduces the width almost with the same accuracy as the TPA, provided that the matching radius $`\overline{r}`$ is large enough to ensure that the contribution from the nuclear attractive potential is small ($`\eta 1`$).

It follows from Table III that $`\eta `$ controls the accuracy of MTPA rather well, except for a broad neutron resonance 1$`f_{5/2}`$ when $`\mathrm{\Delta }\mathrm{\Gamma }/\mathrm{\Gamma }`$ reaches $`10\%`$ at $`\overline{r}=\overline{R}+2.96`$. In this case, the matching radius is quite far away from the barrier radius, so that the accuracy of the MTPA is given by Eq. (11) (with $`R`$ replaced by $`\overline{r}`$). This is well confirmed by Table II, which shows the corresponding correction term.

## V Comparison with the R-matrix theory

It is interesting to compare the final expressions of the MTPA, Eqs. (18) and (23), with the results of the R-matrix theory lane . In the latter method, the space is divided into internal and external regions by a hard sphere of the radius $`a_c`$, and a complete set of the internal wave functions $`u_\lambda (r)`$ is introduced, $`_0^{a_c}u_\lambda (r)u_\lambda ^{}(r)𝑑r=\delta _{\lambda \lambda ^{}}`$. The internal wave functions obey real boundary conditions for the logarithmic derivative:

$$\frac{u_\lambda ^{}(a_c)}{u_\lambda (a_c)}=B.$$
(27)
The value of $`B`$ determines the R-matrix energy eigenvalues $`E_\lambda `$. It is convenient to choose B so that one of the eigenvalues would coincide with the position of the resonance where the value of the phase shift is equal to $`\pi /2`$. This defines the so-called “natural” boundary condition:

$$\frac{u_\lambda ^{}(a_c)}{u_\lambda (a_c)}=\mathrm{}\left[\frac{G_{k_\lambda }^{}(a_c)}{G_{k_\lambda }(a_c)}\right].$$
(28)
In this case, one can apply the one-level approximation, i.e., approximate the resonance with a single Breit-Wigner term. The corresponding width $`\mathrm{\Gamma }_\lambda `$ becomes:

$$\mathrm{\Gamma }_\lambda =\frac{k}{m[1g_\lambda (a_c)]}\frac{u_\lambda ^2(a_c)}{|G_k(a_c)|^2},$$
(29)
where $`g_\lambda `$ denotes the energy derivative of the level shift lane ; erich1 .

One finds that Eqs. (18) and (23) of the MTPA formally resemble Eqs. (28) and (29) by choosing $`a_c=\overline{r}`$ and taking $`g_\lambda (a_c)=0`$. Yet, the inner wave function $`\mathrm{\Phi }_0(r)`$ of the MTPA is different from the internal wave function $`u_\lambda (r)`$ of the R-matrix theory. The latter is totally confined inside the inner region, whereas $`\mathrm{\Phi }_0(r)`$ is a true “bound state” wave function of the inner potential $`U(r)`$. Therefore, their normalizations are different.

The essential problem of the R-matrix theory is a proper choice of the matching radius $`a_c`$, which remains a free parameter. A different choice of $`a_c`$ does affect the results of R-function calculations in a one-level approximation. Moreover, there should exist an optimal matching radius, for which the results of the R-matrix calculations are close to the exact results erich . However, except for some simple cases (e.g., the square well potential), the optimal matching radius cannot be simply prescribed. In contrast, the MTPA is not sensitive to the matching radius $`\overline{r}`$, provided that it is taken inside the “window” defined by Eqs. (11) and (13). This is an essential advantage of MTPA over the R-matrix method. (For critical discussion of the R-matrix expression for the resonance width, see Refs. Arima ; barmore ; kruppaproc .)

## VI Summary

This paper contains a detailed investigation of the two-potential approach to the one-body tunneling problem. It has been found that TPA becomes extremely accurate if the separation radius, dividing the entire space into the inner and the outer regions, is taken deeply inside the barrier, but not too close to the outer classical turning point. From a minimization of the leading correction terms, we obtained simple expressions for the upper and lower limits of the TPA separation radius. The high accuracy of the method was demonstrated explicitly by a detailed comparison with Gamow resonances of a realistic nuclear potential.

Furthermore, we have found that the TPA can be further simplified by taking into account the properties of regular and irregular solutions of the Schrödinger equation under the barrier. The final expressions of the modified two-potential approach formally resemble those of the R-matrix theory with the “natural” boundary conditions. However, the internal wave function of the MTPA is considerably different. In addition, contrary to the R-matrix theory, the corresponding matching radius of the MTPA is well defined. This makes MTPA particularly suitable for practical applications.

## VII Acknowledgments

One of the authors (S.G.) is very grateful to E. Vogt for inspiring discussions and important suggestions. S.G. also thanks TRIUMF, Vancouver and JIHIR for kind hospitality. This work was supported in part by the U.S. Department of Energy under Contract Nos. DE-FG02-96ER40963 (University of Tennessee), DE-AC05-00OR22725 with UT-Battelle, LLC (Oak Ridge National Laboratory), and DE-FG05-87ER40361 (Joint Institute for Heavy Ion Research). T.V. acknowledges support from the Hungarian OTKA Grants No. T029003 and T037991. The support from the University Radioactive Ion Beam consortium (UNIRIB) is gratefully acknowledged.