© ESO, 2014

1. Introduction

In Hillen et al. (2013, hereafter Paper I), we reported the first detection of a large amount (35% of the total) of optical scattered light emerging from the close circumstellar environment of the post-asymptotic giant branch (post-AGB) binary 89 Herculis. Many post-AGB objects display large-scale structures in scattered light (Sánchez Contreras et al. 2010; Siódmiak et al. 2008; Bujarrabal et al. 2002; Ueta et al. 2000), but generally the post-AGB stars in binary systems do not (and are therefore classified as “stellar” in e.g., Siódmiak et al. 2008). These post-AGB binaries (with P ~ 100–3000 d) do have a circumstellar dust reservoir, but of relatively low mass, in the form of a stable Keplerian circumbinary disk (see, e.g., the census of Galactic sources by de Ruyter et al. 2006). This disk is presumably very similar in structure to protoplanetary disks found around young stars, despite the different formation mechanism. In contrast, post-AGB stars that have pronounced aspherical nebulosities, often accompanied by an overdense, dusty equatorial waist or torus that obscures the central object (see, e.g., the Hubble Space Telescope (HST) survey of Sahai et al. 2007), are dubbed protoplanetary nebulae (PPNe). Few binaries, if at all, are found at the heart of PPNe (Hrivnak et al. 2011).

A notable case within the binary group is the Red Rectangle (RR), which has been imaged in various optical filters, showing an extended nebula with a complex spatial as well as chemical structure (Vijh et al. 2006, 2005, 2004; Cohen et al. 2004). The RR is also the only post-AGB binary for which the circumbinary disk has been directly imaged (Lagadec et al. 2011; Cohen et al. 2004; Tuthill et al. 2002; Bregman et al. 1993, and references therein), and for which the Keplerian velocity field has been directly resolved in rotational CO lines (Bujarrabal et al. 2005). The recent Atacama Large Millimeter/submillimeter Array (ALMA) maps of Bujarrabal et al. (2013b) also show a clear molecular component that coincides with the bipolar outflow seen in optical light.

A molecular outflow, in the form of an extended hourglass-like structure, was also detected in the 89 Her system by Bujarrabal et al. (2007). Recently, Bujarrabal et al. (2013a) presented evidence for the presence of low-velocity molecular outflows in a large sample of post-AGB binaries from the shape of the wings in single-dish rotational CO lines. Bipolar outflows might thus be a common property of post-AGB binaries with disks. But it is unlikely that all of these outflows show the complexity of the RR. There is no evidence for a mixed C- and O-based chemistry in most objects as is the case for the RR (Vijh et al. 2006; Waters et al. 1998). The fact that no outflows have yet been detected in optical light in other objects may also be due to a lack of spatial resolution in previous resolved observations (e.g., those of Siódmiak et al. 2008).

Many questions remain about these outflows. How and from where in these complex systems are they driven? There are now detections of jetlike outflows in several post-AGB binaries, thanks to recent and extensive spectroscopic monitoring campaigns (Thomas et al. 2013, 2011; Gorlova et al. 2012; Witt et al. 2009), which suggest an origin in the neighborhood of the companion star. A compact continuum-scattered-light component could go unnoticed in these spatially-integrated spectra, depending on the inclination of the system.

Given the pole-on viewing angle (i ~ 12°), it is possible that our detected 35% optical scattered light around 89 Herculis in Paper I, is largely arising in an outflow. Alternatively, the scattered light could be fully originating from the inner rim of the circumbinary dust disk. The main argument in favor of the latter interpretation is the conformance of the sizes at optical and near-IR wavelengths, as derived from the interferometric data presented in Paper I. The baseline coverage of these data is insufficient to directly distinguish between these two proposed origins: 1) the flux is scattered into our line of sight over an extended surface area of which the outer radius equals the inner rim of the disk; in this case the material is most likely positioned above the orbital plane (i.e., in an outflow or a halo) as the dust would otherwise sublimate; or 2) it arises from a very thin ring, corresponding to the inner rim of the dust disk as resolved in the near-IR. In any case, the observations presented in Paper I require that the projected surface area is more compact at optical than at near-IR wavelengths, so the bulk disk surface can be ruled out as a dominant scattering source.

One way to differentiate between these options is to model the energetics of the disk in detail, taking all spatial constraints into account. If the disk is responsible for all of the circumstellar optical and IR flux, then it must have a significantly larger scale height than if an additional component is present in the system. With an inclination of only ~12°, our line of sight to the disk is moreover unfavorable, given the preference of grains to scatter in the forward direction (see also Mulders et al. 2013). Based on these simple arguments, the presence of an outflow in addition to the disk seems natural. An exploration of full radiative transfer disk models is still needed, however, to rule out the possibility that a disk alone is capable of explaining the data.

The goal of this paper is therefore to investigate whether a self-consistent radiative transfer model of the circumbinary disk in the 89 Herculis system can be created that reproduces all our observations. First, we present some complementary archival mid-IR interferometric and new far-IR photometric observations in Sect. 2. Subsequently, the modeling strategy is described in Sect. 3, and a first confrontation with the observations is done in Sect. 4. In Sects. 5–7, some improvements are made to the model, and the fit to the near-IR, mid-IR and optical interferometric data is examined in detail. Finally, in Sects. 8 and 9 we discuss our results and summarize our conclusions.

2. Observations

A detailed account of our optical and near-IR interferometric data was given in Paper I. Here, we add some archival data from the MID-infrared Interferometric instrument (MIDI) on the Very Large Telescope Interferometer (VLTI) to the existing data set, as well as new Spectral and Photometric Imaging Receiver instrument (SPIRE) photometry from the Herschel satellite.

2.1. Additional interferometry

Three VLTI/MIDI observations of 89 Her, procured in April 2006 and already discussed in Bujarrabal et al. (2007), were reduced with version 2.0 of the EWS software (Jaffe 2004)¹. This software estimates the MIDI visibilities in a coherent and linear way, based on a shift-and-add algorithm applied on individual interferometric frames. Two observations used Auxiliary Telescope (AT) stations D0 and G0, and the other observation used stations E0 and G0. The projected baseline angles are 75–90° east of north, and the lengths are ~30 and 15 m.

In the calibration, the raw target visibility V_sci is divided by the instrumental visibility V_ins for each wavelength: (1)Here, C and P are the measured correlated flux and photometric counts for the target (sci) and calibrator (cal), and V_{theo − cal} is the expected calibrator visibility. This assumes that the calibrator and target observations are performed under the same instrumental conditions. For a robust calibration, we interpolate the V_ins of the calibrator directly preceding and succeeding the target observation. To account for the longer term variability of the transfer function, we add the standard deviation of V_ins for all calibrators within four hours of the science measurement in quadrature to the visibility errors. On the worst night, the variability of the transfer function is ~20%, which can (partly) be attributed to the difficulty in estimating the photometric counts, given the strong background radiation of the delay-line system. This noise source might be source-dependent because of the adjustments in the delay lines upon a change in the coordinates on the sky.

The calibrated visibilities are shown in Fig. 1. The 30 m baseline observations are partly resolved at the blue edge of the band, with visibilities of V ~ 0.8. At the red edge, the situation is less clear due to the larger errors, but both 30 m observations seem to be in agreement with a slightly resolved object. The single 15 m baseline observation shows a flat visibility curve with V ~ 0.8. This suggests that 89 Her is at least as resolved on the 15 m baseline as on the 30 m baseline, however, we are cautious in our interpretation because the observing conditions were poorer during the night of the 15 m observation than for the other two nights (with an average atmospheric coherence time of 2.5 rather than 6 and 11 ms). The resulting stronger optical path difference fluctuations make it more difficult to avoid decorrelation effects when reducing the data², which may effectively lower the apparent correlated flux. Given the stability of the transfer function during the night, we do not exclude the 15-m baseline, but keep in mind that the observation is potentially affected by the poorer atmospheric conditions.

2.2. PTI fluxes

In addition to visibilities, fluxes can be extracted from the raw data of the Palomar Testbed Interferometer (PTI, Colavita et al. 1999). We developed a calibration procedure for these fluxes, which has been detailed in Hillen et al. (2012). The PTI measured visibilities in five channels within the K band and four channels within the H band. Only the central channels within the H and K bands are used in our interferometric analysis, but we include all channels for the flux calibration. As explained in Hillen et al. (2012), we separate the calibration into an absolute calibration in one of the channels within each band (at 1.7 μm and at 2.2 μm, where the atmospheric transmission is the highest) and a separate calibration of the spectral shape over the other channels in the bands. This nonstandard approach was developed since the spectral shape calibration can be done more accurately, even in less optimal atmospheric conditions, and lead to accurate light curves of Mira variables (see Hillen et al. 2012). For 89 Herculis the calibrated fluxes are also stable in time and are averaged to obtain the final fluxes that are shown in Fig. 2. These data will serve as an independent check of the shape of the spectral energy distribution (SED) in the near-IR spectral range and how the best models perform at these wavelengths.

Fig. 1 Three baselines of MIDI data overplotted (at a single continuum wavelength, indicated in cyan) with all double power-law MCMax models computed in Sect. 5. The dotted blue line, “Best 1” in Table 1, is the best single power-law MCMax model in terms of fitting the near-IR visibilities (see Figs. 4 and 5). The dashed red line corresponds to the best-fitting double power-law model with the metallic Fe opacity source (“Best 2” in Table 1). The full green line shows the best-fitting double power-law model when metallic Fe is replaced by amorphous C (“Best alt.” in Table 1). See Sects. 4.2, 5.1, and 6 for a discussion of the models. The inset shows the data at the 30 m baselines as a function of wavelength, with the same double power-law models overplotted as in the main panel.

2.3. SPIRE photometry

The object 89 Her was observed with the SPIRE instrument (for more info on its scientific capabilities, observing modes, data reduction, calibration, and performance, see Griffin et al. 2010; Swinyard et al. 2010) aboard the Herschel satellite (Pilbratt et al. 2010) in March 2012 with a 169 s exposure time (Herschel Observation ID 1342240026). The observation was performed simultaneously in three passbands (centered at 250 μm, 350 μm, and 500 μm), covering an area of 8′ × 8′ homogeneously. The naiveMapper task was used to project each observed sample onto a single pixel in the reconstructed map. The pixel size is 4′′, 6′′, and 9′′ in the three bands, respectively, and the corresponding beam FWHM are 18.2′′, 24.9′′, and 36.3′′, respectively. Aperture photometry was used to extract the final fluxes, after sky background subtraction. The latter was estimated with the standard daophot algorithm. The errors are dominated by the absolute flux calibration and were assumed to be 15% (Swinyard et al. 2010). The final calibrated fluxes are shown with the rest of the SED, in Fig. 3.

2.4. Data selection

Our photometric data were already discussed in Paper I, with the exception of the Plateau de Bure Interferometer (PdBI) 1.3 and 3 mm continuum fluxes, which we have taken from Bujarrabal et al. (2007). We decided to exclude the Infra-Red Astronomy Satellite fluxes (IRAS, Neugebauer et al. 1984) from our χ² calculation, as we found that they are systematically higher than the AKARI (Murakami et al. 2007), Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010), SPIRE and Submillimetre Common-User Bolometer Array (SCUBA Holland et al. 1999) measurements at similar or longer wavelengths. Since the Infrared Space Observabory (ISO) spectrum is very noisy beyond 25 μm, as well as to simplify the analysis, we decided to only include it by integrating it over several well-chosen photometric passbands in the near- to mid-IR domain. These passbands are chosen so as to fill up the voids in the spectral coverage of the real photometry, and are the Infrared Array Camera bands (IRAC, aboard the Spitzer Space Telescope, Reach et al. 2005) at 5.8 and 8.0 μm, as well as the WISE W2 (~5 μm) and W3 (~12 μm) bands, and the Midcourse Space Experiment (MSX, Mill et al. 1994) D band at 15 μm. Since the error on the absolute calibration of the ISO spectrum is of the order of 10%, this error value is also adopted for the photometry derived from it.

Fig. 2 PTI absolute fluxes (in the insets) and spectral shapes over the H (upper panel) and K (lower panel) bands. The data are shown in blue. The dashed red and full green lines correspond to the same models as in Fig. 1.

Fig. 3 Spectral energy distribution of 89 Her with the same three models as in Figs. 1 and 2 overplotted. The dot-dashed gray line shows the input Kurucz model, here reddened with an interstellar contribution of E(B − V) = 0.07 (see Paper I). The black and gray data points are the “stellar” V,R,H, and K band fluxes as determined in Paper I. Open symbols represent actual measured fluxes in the passbands that are listed in the legend. The full black line is the ISO spectrum. The dotted blue line, “Best 1” in Table 1, clearly does not agree with the observed SED. A better, but still not excellent, fit to the SED can be obtained with an inner radius twice as large, in combination with p = 1.5, but this results in first-lobe near-IR visibilities that are too small (see the text). The dashed red and full green line correspond to the best-fitting double power-law models with metallic Fe and amorphous C opacity sources, respectively (see Fig. 1). The MCMax models have been reddened similarly as the Kurucz model.

3. MCMax modeling

3.1. Introducing MCMax

The 2D Monte Carlo radiative-transfer code MCMax was used to solve the structure of the disk and subsequently calculate model SEDs and visibilities. The MCMax code has been successfully used to model both spatially and spectrally resolved observations of different kinds of gas-rich disks (Mulders et al. 2011; Verhoeff et al. 2011; Acke et al. 2013). Our modeling approach was basically the same as presented in Mulders & Dominik (2012) and we refer to this paper and Min et al. (2009) for a full account of the procedures and physics implemented in MCMax. In short, MCMax combines the immediate re-emission procedure of Bjorkman & Wood (2001) with the method of continuous absorption by Lucy (1999) to compute the radiative equilibrium temperature stratification throughout the disk. Subsequently, the vertical scale height of the disk is self-consistently determined by solving the equations of vertical hydrostatic equilibrium, implicitly assuming thermal coupling between the gas and dust, but dynamical decoupling, i.e., dust settling. The final disk structure is found by iterating both steps until convergence is reached, typically in less than five iterations. Scattering by dust can be treated with a full angle-dependent Mueller matrix for each scattering event, or by assuming isotropic scattering, which needs less computation time (Min et al. 2009).

All observables are calculated by analytically integrating the computed source function by ray-tracing. This includes the visibility: the Fourier transform is analytically calculated on a circularly sampled image. The advantage of this approach is that the inner parts of the image are sampled at a higher resolution so that visibilities can be reliably computed up to high spatial frequencies.

3.2. The model setup

We define a “standard” radiative-transfer disk model based on some general characteristics for the sample of disks around post-AGB binaries:

• Deroo et al. (2007) made a diskmodel of IRAS08544-4431 by fitting the total disk mass and innerrim radius to the SED, while fixing all other parameters to typicalvalues for the group of binary post-AGB objects or based on otherexternal constraints (e.g., thestellar T_eff from optical spectroscopy). By tuning the inclination and position angle, this “standard disk model” reproduced all of their interferometric observables surprisingly well.

• All disks around post-AGB binaries (with a few clear exceptions) have an O-rich mineralogy, their mid-IR spectral features are dominated by amorphous and crystalline silicates and their dust typically has a high degree of grain processing (Gielen et al. 2008, 2011).

• All objects for which long-wavelength fluxes are available in the sample study of de Ruyter et al. (2006) have blackbody slopes that are indicative of the presence of a component of mm-sized grains.

• Sahai et al. (2011) presents evidence for the presence of grains >1 mm in the compact disks of RV Tau, U Mon, and AC Her, from flux measurements at sub-mm and mm wavelengths. The derived dust masses lie between 5 × 10^-4 and 1 × 10^-3 M_⊙.

The overall shape of the SED of 89 Her is very reminiscent (see also de Ruyter et al. 2006) of a group 2 Herbig star as defined by Meeus et al. (2001). We therefore include dust settling in our disk model in a self-consistent manner (see Mulders & Dominik 2012, for more details about the implementation and consequences of this approach). Our dust composition and grain shape are the same as in the latter paper (i.e., irregular shapes, approximated by the distribution of hollow spheres (DHS), Min et al. 2005), except that we start with metallic Fe as the additional featureless “continuum” opacity source (~15% in mass fraction) instead of amorphous carbon, which would be more difficult to reconcile with the oxygen-rich circumbinary environment. The majority of the dust mass consists of silicates (for the composition see Min et al. 2007).

The dust is distributed over a range of grain sizes with a power-law dependence, f(a) ∝ a^{− q}, with the power-law index q as a fitting parameter. We attempted the values 2.5, 3.0, 3.25, and 3.5. The minimum and maximum grain sizes are set to 0.01 and 10⁴ μm. The smallest grains are needed to have sufficient opacity at the stellar maximum, while the inclusion of large grains is based on the findings listed above, and the confirmation that for 89 Her there is also a non-negligible mm flux. To implement the grain-size-dependent settling in practice, the continuous size distribution is divided into a finite number of bins. The grain radius of each bin is defined as the logarithmic mean of the minimum and maximum grain size within the bin. Two bins are used per order of magnitude in grain radius, and in each bin the opacities are calculated as the average over 10 sizes to avoid introducing strong resonant features. The turbulent mixing strength is fixed at the standard value of 0.01 (Mulders & Dominik 2012). The same goes for the gas/dust ratio, which is taken to be 100.

One of the largest unknowns in modeling disks is the radial dependence of the surface density distribution, which is typically parameterized with a simple power law Σ ∝ r^−p. Initially, we do not deviate from the single power-law formalism, but do explore a range of values for the index p. Unlike for protoplanetary disks where typically p is found to be ~1.0 (see Williams & Cieza 2011, for a general overview), in the post-AGB case this index has not yet been constrained with resolved millimeter observations. Moreover, there is no a priori reason to assume that it is the same as in protoplanetary disks, given the difference in formation scenario.

The inner and outer radius of the disk, as well as the total dust mass are also free parameters of our model. For the inner radius, small step sizes of 1 AU are used because our interferometric data are sensitive to this parameter. A firm upper limit on the disk’s outer radius comes from the limiting resolution of 400 mas in the PdBI maps of Bujarrabal et al. (2007). The diffraction-limited images of Lagadec et al. (2011) impose a thermal IR upper size limit of 300 mas. The HST scattered light images of Siódmiak et al. (2008) have a resolution of ~75 mas. The two latter constraints do not imply that there is no material beyond the respective radii, it simply means that any material that might be there only emits/scatters with a very low contrast at the given wavelengths. Nevertheless, these results do indicate that the disk is likely compact. The effect of the outer disk radius on the SED is rather small, and we only vary it because it affects the surface density distribution in a different manner than changing the power-law index.

Finally, the stellar properties are also required for the radiative transfer calculations and to determine the vertical structure of the disk. Some of these values were discussed extensively and redefined in Paper I. The stellar mass is fixed to 1  M_⊙, taken as the sum of a 0.6  M_⊙ post-AGB star and a 0.4  M_⊙ companion, so that the corresponding inclination is ~13.5° (derived from the measured mass function of Waters et al. 1993). To correct for interstellar reddening, we use the reddening law of Fitzpatrick (2004) and an E(B − V) = 0.07 (see Paper I). The final list of input parameters and associated values are listed in Table 1.

4. First results

In this section, we perform a first confrontation between the radiative transfer models and the observations, from which it will become clear that our model is too simple to fit the observations. In subsequent sections attempts will be made to resolve the discrepancies.

4.1. An initial selection

A large grid of MCMax models was calculated with the input parameters listed in Table 1. Our fitting strategy consisted of several steps. First, for computational efficiency, the full grid was calculated under the assumption of isotropic scattering, and without a full ray tracing of the spectrum. In the case of isotropic scattering, the raw Monte Carlo spectra can be used as a good proxy for the real spectrum, as long as sufficient photons are used. For comparison with the photometric observations, we integrate the model spectra to compute synthetic photometry for each of the passbands in our data set, taking into account whether a bolometer- or CCD-type of detector is involved.

Each model in the grid is then first compared to the IR part of the SED. To this end, we compute a normalized χ² that includes all photometry at wavelengths larger than 6 μm. Only mid- to- far IR fluxes are included in the χ² because at these wavelengths the effect of scattering is small. In this way, a first selection can be made of models that generally fit the observed IR-excess. A very rough cutoff criterion of was applied. The advantage of this approach is that the correlations that exist between different parameters, and that cannot be resolved by SED-fitting, are withheld. From inspection of the retained models, we note:

• Most models that fit the χ²-criterion have a dust size distribution power-law index q of 3.0 or 3.25, but a minority of models with 2.5 and 3.5 cannot just be excluded.

• The q = 3.5 models generally underpredict the (sub-)mm fluxes, tend to give a too strong 10 μm feature, and have too large near-IR fluxes.

• For q = 2.5, the far-IR slopes are too shallow, and dust masses larger than 10^-3 M_⊙ are required.

• The best-fit models with q = 3.0 or q = 3.25 have a dust mass in the range 10^-4–10^-3 M_⊙.

• The outer radius of the disk is not constrained, although models with radii of 50 and 100 AU more often lead to good χ² values.

• The inner radius and the power-law index of the surface density distribution are not constrained individually, as expected.

Bujarrabal et al. (2007) derived a total gas mass of ~10^-2 M_⊙ for the central (disk) component that was unresolved in their CO-maps, assuming a distance of 1 kpc. This was recently confirmed from single-dish CO observations by Bujarrabal et al. (2013a). Since we assume a larger distance (see Paper I), and because of the inherent uncertainties in estimating the total gas content from ¹³CO lines, we allow our models to have a maximal total mass of 0.1 M_⊙, or equivalently a dust mass of 10^-3 M_⊙. The minimum dust mass needed to fit the SED is 10^-4 M_⊙.

4.2. A full scattering treatment

All models that fit the initial χ² criterion, and the imposed mass constraints, are now recomputed with full scattering. A ray-traced spectrum is now needed, as in the full angle-dependent scattering case, the Rayleigh-Jeans tail of the Monte Carlo spectra becomes unreliable. To compute the visibilities, a disk position angle of 0° east of north is assumed. Given the small inclination, this value cannot be determined as it has a very small influence. A good check on the amplitude of the potential error, is to consider the effective baseline defined by Tannirkulam et al. (2008): (2)with B_proj the actual projected baseline length and θ the angle between the projected baseline and the major axis of the disk. With our assumed inclination, the correction term amounts to at most 2%.

Different χ² are defined as a measure for goodness-of-fit. The SED is split into three wavelength domains, based on their dominant circumstellar emission mechanism: 1) scattering dominated for wavelengths below 1.5 μm; 2) thermal emission dominated beyond 6 μm; and 3) a near-IR transition regime where both mechanisms are equally important and where our visibilities impose strong spatial constraints. A different χ² is defined for each of these domains. Similarly, separate χ² values are computed per passband (R,H,K, and N) with measured visibilities. The performance of the model at optical wavelengths will be evaluated in Sect. 7. First, we search for models that fit all IR constraints, i.e., the IR-excess as well as the mid- and near-IR visibilities.

Fig. 4 H band interferometric data. The left and right panels show the H band visibility amplitudes and squared visibilities observed by the various instruments. The visibility amplitudes, as produced by the data reduction software packages of the CLIMB and CLASSIC instruments, are not squared in order to preserve the original error statistics. The data are the same as in Paper I. The same three models are shown as in Figs. 1 and 3 and as listed in Table 1.

Fig. 5 Same as Fig. 4 but for the K band.

We select good models based on the following criterion: for all of the data sets. Taking a minimal of 1.0, allows us to avoid excluding too many models in case a particular data set is over-fitted. On the other hand, not imposing throughout allows us to find models that globally fit all data sets, without necessarily fitting each data set to high precision. One must keep in mind that there is some variability in the data that cannot be reproduced by our model. In particular, the near-IR visibilities show some scatter that is larger than the errors of some of the corresponding data, which could be a time-dependent effect (perhaps due to the binarity), the result of spatial/azimuthal flux variations (e.g., a low-amplitude spiral pattern in the disk) or just an underestimation of some of the errors.

None of the current models fits our criterion. A closer inspection of the best-fitting models shows that the near-IR visibilities cannot be fitted simultaneously with the SED. Basically, the models that generally fit the SED well have an inner rim radius that is too large to fit the first lobe of the H and K visibility curves (which dominate the χ²). In turn, the very few models that fit these first lobes better, and therefore come out with good visibility χ²’s, overpredict the H band flux while underpredicting the near- to mid-IR fluxes, and therefore do not fit the intermediate nor the long near-IR baselines either (see the blue dotted line in Figs. 3–5; and “Best1” in Table 1). Models with small inner radii (~3.5–4.5 AU) that do predict enough flux beyond ~K band were already removed in the previous section, as these usually had extreme near- to mid-IR excesses when computed with isotropic scattering. As such, the remaining models with small inner radii generally are those with flat surface density distributions. As a check of our strategy, we recompute the 20 best models with an inner radius of 4.5 AU that were previously removed with a full scattering treatment, and find that these models do not fit our criterion.

Fig. 6 Mass absorption (left) and total extinction (right) coefficients of the mixtures used in this work. Our “standard” mixture (full lines) contains 15% metallic iron, while the alternative mixture (dashed lines) has 15% amorphous carbon instead. The latter composition is exactly the same as that used by Mulders & Dominik (2012). Only six grain sizes are shown for clarity.

5. A good fit to the IR data

To improve on the structural deficiency of our model, we change our surface density distribution parameterization and adopt a double power-law formalism (see Sect. 8 for a further motivation and discussion of this choice). Since this new parameterization introduces two more fitting parameters, the turnover radius R_mid and the inner density power-law index p_in, we reduce the parameter space for four other parameters: the outer radius and surface density power-law index (p_out), the total dust mass, and the grain size power-law index q. Given the indications for a very compact disk with large grains (see Sect. 3.2), but only a moderate mass (see Sect. 4.1), and since a significant scale height is required to explain our optical observations, we adopt our most compact configuration of R_out = 50 AU. Furthermore, we adopt the “canonical” value of p_out = 1.0 or 1.5. The dust mass is limited to 1.0, 2.5 or 5.0 × 10^-4 M_⊙. The value of q can be either 3.0 or 3.25. Based on our geometric analysis of Paper I, and the good correspondence with the first-lobe visibilities of models with inner radii of ~4–5 AU in the previous section (see Figs. 4 and 5), we vary the inner radius between 3.75 and 4.75 AU, in steps of 0.25 AU.

5.1. The SED and near-IR visibilities

The good models are searched in the same way as described in the previous section. Five models in the grid now conform to our acceptance criterion, and in general the χ²’s have decreased significantly below the values of the single power-law models. Interestingly, all good models have q = 3.0 (imposed by the H band visibilities) and a dust mass equal to 5 × 10^-4 M_⊙. There are no clear preferences for the surface density gradients. The inner and turnover radius are also not well defined, but this can be expected since we already refined the sampling in these parameters. In any case, the turnover radius lies about three times further than the inner radius of the disk on average.

In the 2–6 μm range, a flux deficit is present when compared to the ISO spectrum, however, based on the photometry, the fit and the corresponding χ² is very good. The H and K band fluxes fit quite well, as is confirmed by the PTI fluxes in Fig. 2. The slope of the spectrum does not exactly agree: there is slightly too much flux in the long-wavelength end of the H band and too little flux in the entire K band. Also, the spectral slope towards 2.4 μm is too shallow with respect to the data, which confirms the mismatch with respect to the short-wavelength end of the ISO spectrum. We come back to these small discrepancies in Sects. 6 and 8.

The model fits the visibilities (shown in Figs. 4 and 5) quite well in both the H and K band. The rather high χ² in the H band suggests otherwise, but this can largely be attributed to the previously mentioned scatter in the data. The largest contribution to the H band residuals comes from the IOTA data around spatial frequencies of 115 cycles/arcsec and the PTI data near 325 cycles/arcsec.

The parameter values of the preferred model are listed in Table 1 as “Best 2,” and we accept it as the model that fits our IR data sets best. The structure of this model will be examined in more detail in Sect. 8.

5.2. The mid-IR visibilities

The MIDI data have an insufficient spatial resolution to strongly constrain the position and shape of the inner disk rim, and the data are more sensitive to the bulk of the inner disk. Therefore, it is surprising that the visibilities (see Fig. 1) are so similar at the 15, and 30 m baselines. The predictions of all our double power-law models and of the “best” single power-law model are depicted as well. Clearly, none of the models simultaneously fits both baselines, or even just the slope between them. Nevertheless, basically all of the models fulfill our acceptance criterion and have a good χ², thanks to the larger weight of the 30 m baseline (as it is observed twice) and the rather large errors of the data. The shape of the MIDI visibility spectrum does not hold much information, given that this shape varies between the two 30 m observations. While our best-fit double power-law MCMax model is not perfect, it formally reproduces the mid-IR data sufficiently. In Sect. 8.3, we return to the 15 m discrepancy.

6. Amorphous carbon?

Finally, we also computed a grid of models, with the same parameter values as in the previous section, but with 15% of amorphous carbon instead of metallic Fe as the additional opacity source (see Sect. 8.4 for a discussion of our motivations). As is illustrated in the left panel of Fig. 6, amorphous C has a higher mass absorption coefficient in the near- to mid-IR domain compared to metallic Fe, independent of the grain size, and in the (sub-)mm domain for the large grains. In terms of total extinction (as shown in the right panel), however, the differences are much smaller, except for the smallest grain sizes. Therefore, one might expect a better agreement of these models to the ISO spectrum and the exact shape of the SED in the near-IR.

We applied the same strategy to find the best-fit models within this grid (see also Table 1 and Figs. 1–5). The near-IR fluxes are reproduced marginally better, as are the mid-IR fluxes, the silicate feature at 10 μm, and the H band visibilities. On the other hand, the fit to the (sub-)mm fluxes degrades. The reason for these changes lies in the changed values of several parameters: the models are now allowed to have relatively more small grains (q = 3.25), but only if the total dust mass (M_dust = 1 × 10^-4) decreases accordingly. The surface density is now also generally flatter with p_out = 1.0.

Fig. 7 NPOI and VEGA squared visibilities at 673 nm. Our complete grid of double-power-law MCMax models is also displayed (see Sect. 5). The best-fit models listed in Table 1 are shown with the same color coding as in Figs. 1–5.

7. Fit to the optical data

We finally confront the models with the optical visibilities, and the optical part of the SED. Fig. 3 already illustrated that the best IR models of the previous sections fail at wavelengths below 1 μm. This is the case for all of the calculated models. Of all the models computed in Sect. 5, the model that performs best still has a χ² ~ 80 for the optical part of the SED, and underpredicts the observed scattered flux. To illustrate this more clearly, Fig. 7 shows the NPOI and VEGA visibilities at 673 nm, with all of the full-scattering models that were calculated in Sect. 5 overplotted.

Not a single model reaches the low visibilities that were observed. Interestingly, the best model of Sect. 4.2 (the blue dotted curve on Figs. 3 and 7) performs better than our best double power-law models computed in the previous sections. Models with amorphous C generally perform worse at optical wavelengths than models with metallic Fe, due to their smaller single-scattering albedo (see Fig. 8) and larger fraction of small grains. For our metallic Fe models, much higher fractions of optical scattered flux are found on average with single rather than with double surface density power laws. These models have large dust masses, however, often in combination with a large fraction of small grains (q = 3.25), so they are totally incompatible with the IR constraints as they predict too much H, K and often even mid-IR, fluxes. The best double power-law models, with a VIS_optχ² of only 3.1, suffer in fact from the same problem, which is most pronounced in the H band. A closer inspection shows that the lowest χ² actually does not correspond with the highest amount of optical scattered light, which is reserved for the models with large turnover radii R_mid. The latter however underpredict the visibilities at the short spatial frequencies, which have a larger weight in the χ² simply by their number, in particular in the channels near 0.85 μm (not shown in Fig. 7). Basically, these models have a larger area over which the disk surface contributes to the scattered light (see also the τ_V = 1 surface in Figs. 9 and 10). Since the short spatial frequencies are fairly well reproduced by the majority of the models, the resulting range in χ² is rather modest: the worst model only has a χ² ~ 22. Finally, we note that many models, including our general best-fit model, result in very flat visibility plateaus at the high spatial frequencies, as is suggested by the data as well.

Fig. 8 Albedo of the smallest three particle sizes and for four different compositions: silicates + amorphous carbon (dashed lines), only silicates (dotted lines), silicates + 15% iron (full lines), and silicates + 30% iron (dot-dashed lines).

Fig. 9 Temperature structure of our best-fit models. The upper panels show maps of the temperature, with the radial and vertical τ = 1 surfaces at 0.55 μm overlaid as solid and dashed lines, respectively. The middle and lower panels display the midplane pressure scale height and integrated dust density, respectively, as a function of the disk’s radius. On the left and right, the best Fe and amorphous C containing models are shown, respectively.

Fig. 10 Density structure of our best-fit models. From top to bottom, maps of the gas density, the dust density, the dust density for particles of ~0.06 μm, and the dust density for particles of ~5 mm are shown. On all panels the radial and vertical τ = 1 surfaces at 0.55 μm are overlaid as well, with solid and dashed lines, respectively. The left and right panels contain the best Fe and amorphous C containing models, respectively.

8. Discussion

In this work, we have shown that with reasonable input parameters, and using a self-consistent treatment of dust settling, radiative transfer disk models can be constructed that are in good agreement with all our IR observations within their uncertainties. The models successfully reproduce the general trend in the near-IR visibilities and fit the full SED up to sub-mm wavelengths very well. On the other hand, these best models, and in fact every one of our models, are incapable of producing a sufficient amount of optical scattered light. Before we discuss the deficiencies, we put our results in perspective.

8.1. The inner rim curvature

That a simultaneous fit to the SED and the near-IR visibilities could not be found with a single power-law model is not surprising. Our near-IR interferometric observations clearly prefer an extended and smooth emission region (see also Paper I), not the wall-like emission profile resulting from a single power-law parameterization of the surface density. For protoplanetary disks, interferometric investigations have already extensively shown that the inner disk rim is usually not sharp, but has a smooth geometry (see Dullemond & Monnier 2010, for a nice review). Different physical processes have been identified to cause this rounding off.

Isella & Natta (2005) showed the importance of including an accurate physical description of dust sublimation and evaporation, and in particular its temperature dependence on the gas pressure. Tannirkulam et al. (2007) for the first time presented a 2D radiative transfer model of a rounded-off rim by considering the effect of dust growth and settling. Due to these effects, big grains, which survive closer to the star thanks to their ability to cool more efficiently, are also primarily found in the midplane, and thereby round off the rim.

The results of Isella & Natta (2005) were later confirmed by a more extensive exploration of the effect of global disk properties on inner rim models, using a detailed physical treatment of dust formation and destruction in 2D radiative transfer models (Kama et al. 2009). An important assumption in the latter work is their double power-law radial parameterization of the surface density, as assumed here. Their smooth inner structures, often having wide optically thin dust regions in front of the actual rim, can partly be attributed to this parameterization, which is probably more realistic than the wall resulting from a single power-law model, as we started with in Sect. 4.

Our models are not computed with a sophisticated dust sublimation/evaporation treatment, like in Kama et al. (2009), since it requires an immense computation time to do so for a grid of models. The self-consistent inclusion of a grain size distribution, combined with dust settling, already increased the computation time of a single model from tens of minutes to several hours (on a single CPU). Kama et al. (2009) found smaller inner rim radii than Isella & Natta (2005) and Tannirkulam et al. (2007). This was attributed to the implementation of an abundance correction to the used dust sublimation law by Kama et al. (2009). It would be interesting to explore how these detailed inner rim curvature calculations exactly depend on the abundances, and what would be the effect in a post-AGB context. In particular, it would be interesting to investigate how our best-fit models would evolve if dust destruction were properly included in the calculations, but this is beyond the scope of the current paper. We also note that Kama et al. (2009) did not use a grain size distribution for their dust, and hence neglected the effect of settling.

Instead of extending the model physics, we preferred in this paper to do a wider exploration of the model parameter space, since post-AGB disks have not yet been studied extensively. In the following subsection we check, a posteriori, the temperature and density structure of our best-fit models.

8.2. The disk structure

Figure 9 shows the temperature structure of our best-fit Fe (left) and amorphous C (right) containing models as a function of the radial and vertical position in the disk, overplotted with the radial and vertical τ = 1 surfaces at 0.55 μm. The figure also shows the radial profile of the gas pressure scale height in the midplane (middle panels) and the radial profile of the vertically integrated surface density (lower panels).

An inspection of the inner rim dust temperatures points out that they are rather high, up to ~1600 K. The model with amorphous C has a slightly smaller temperature at the inner rim, but only by ~50 K. Interestingly, the model with metallic Fe also has a larger radial temperature gradient compared to the amorphous C model. High dust temperatures have appeared in the literature before. Benisty et al. (2010) even needed values up to 2200 K, but these values may not be realistic. Considering the results of Kama et al. (2009), who attempted to implement the most accurate treatment of dust sublimation possible, and adopting our gas density of ≤10^-11.5 g/cm³ (see Fig. 10), sublimation temperatures of only 1300 K for olivine and 1400 K for iron are to be expected. It can therefore be foreseen that by including dust sublimation properly, our inner rim radius would increase, which would be incompatible with the near-IR interferometric data.

The gas as well as dust density structure of the same two models are depicted in Fig. 10. The figure also separately shows the density of small (i.e., 0.06 μm) as well as large (i.e., 5.6 μm) grains, to illustrate the relevance of dust settling. The Fe containing model in general has higher densities than the C models, except in the very small grains. This is expected, given the differences in total dust mass and grain size distribution (see Table 1).

We fixed some parameters and made some implicit assumptions in our modeling. These specific choices could impact our results. Our total dust mass agrees with the value determined by Bujarrabal et al. (2013a) based on single-dish, rotationally-excited CO line fluxes, in particular for our C-containing model. If the true disk outer radius would be larger than the 50 AU assumed in our double power-law models, however, the required dust mass may increase as well. On the other hand, the mass derived by Bujarrabal et al. (2013a) needs to be validated with full radiative transfer calculations of the line profiles, taking the disk structure into account. Another parameter could be the turbulent mixing strength α_turb, which determines the efficiency of vertical mixing, and hence strongly impacts on the vertical structure of the disk. Similarly, the central star’s mass and the global dust/gas ratio influence the disk’s scale height. We calculated a few models with α_turb increased or decreased by an order of magnitude or with M_⋆ = 0.75 M_⊙, and it is unlikely that these parameters will affect our main results.

8.3. The MIDI 15 m baseline

The discrepancy between our models and the 15 m MIDI baseline could be because of the poorer atmospheric conditions during the night of that measurement. If this is not the case, however, then the smoothness of all of our model visibility curves shows that to reproduce both observations, a simple continously declining radial density profile does not work. Instead, a strong decrease in the optical depth of the disk just behind the inner rim is needed. This means that to fit these data, a gap must be developing in the disk between ~10 and 40 mas in radius, or that there must be a strong shadow that is not present in our MCMax models.

There is an ongoing debate about shadows and puffed-up inner rims in protoplanetary disks (see, e.g., the review of Dullemond & Monnier 2010). In some of the original models (e.g., those of Dullemond & Dominik 2004), the region immediately behind the inner rim had a smaller scale height than expected from a fully flaring disk. The interpretation was that this was because of a strong puffing-up of the frontally illuminated inner rim, which then casts a shadow on the region immediately behind it. It was, however, questioned whether such a shadow could lead to a (partial) collapse of the outer disk (see, e.g., the discussion in Wood 2008). Instead, it was suggested that the very steep surface density gradient in the self-shadowed models is responsible for the disk not reaching the large scale heights of flaring disks. Also, the too simple treatment of scattering and radial radiative diffusion of heat might have played a role in the early models, which is no longer a problem in the more recent Monte Carlo based radiative transfer codes.

Based on the Dullemond & Dominik (2004) models, van Boekel et al. (2005) presented simulated visibility curves at 10 μm of a flaring (with partial “shadowing”) and fully self-shadowed case. Interestingly, these curves have a plateau in the continuum at 7.7 μm because of the presence of a shadow in the models. Despite that our observations look very similar in the continuum, the inadequate treatment of scattering and diffusion in these models cannot be neglected. Therefore, if real, a gap interpretation to explain our observations would seem more likely.

The SED does not show any sign of a gap in the disk, while for protoplanetary (pre-)transition disks, in which gaps have been spatially resolved, this is mostly clearly the case (Williams & Cieza 2011). For this reason, we assume that our 15 m observation is untrustworthy, but more mid-IR interferometric observations are needed to confirm this.

8.4. The dust composition

Our choice to explore models with amorphous C was inspired by the recent result of Acke et al. (2013), and by the need to investigate whether a change in dust composition can significantly affect the predicted amount of scattered light by the disk (see Sect. 8.5).

Acke et al. (2013) claim that the continuum opacity source in the circumbinary disk of HR4049 can be identified as amorphous carbon, by confronting radiative transfer models with near- and mid-IR interferometric data and the IR ISO spectrum. They state that the identification can be made in essence because the temperatures reached by various dust species at a given physical distance to the star are very different. As an explanation for the origin of the carbonaceous material, they propose that ongoing mass loss from the depleted photosphere of the post-AGB primary consists of carbon-bearing molecules. This is consistent with various observations (Geballe et al. 1989; Roberts et al. 2012; Acke et al. 2013; Malek & Cami 2014). According to Acke et al. (2013), this material may then form a veil of carbon-rich dust over the pre-existing silicate-rich disk. It remains to be demonstrated, however, that this outflowing material influences the structure and appearance of the disk, such that this is detectable.

The models used by Acke et al. (2013) were rather simple and broke down beyond 10 μm, as the dust consisted of only 0.01 μm grains of a single species. Although the effect of very small metallic Fe vs. amorphous C grains is very different in terms of heating, the influence of the larger grains, which constitute the bulk of the dust mass, on the temperature structure cannot be neglected (see also Mulders & Dominik 2012). Our results show that their conclusion about the composition of the dust component in the HR4049 disk needs to be carefully reevaluated. Our best Fe and C models give very similar fits to our observations of 89 Her; and the latter constraints are similar to those relied upon by Acke et al. (2013). As Figs. 9 and 10 show, the inner temperature structures are not very different between our two best models, while the density structures are different. We stress that the disk structure needs to be taken into account, including the important effect of dust settling, like was done in this work.

No evidence exists for the presence of carbon-rich gas nor dust in the 89 Herculis system. A large-scale outflow has been detected in CO (Bujarrabal et al. 2007), but its origin within the central object is unknown (see Sect. 8.6). Also, the photospheric depletion pattern in 89 Her is not very pronounced (Kipper 2011), in contrast to HR4049, which is one of the most depleted objects (Waelkens et al. 1991) among the known sample of Galactic post-AGB binaries. Hence, from the perspective of abundances, the photosphere of 89 Her contains sufficient Si and Fe to form O-bearing dust in the outflow. While our best-fitting C model matches the observations mildly better than our best-fitting Fe model, it seems much harder to explain the formation of amorphous carbon dust in the 89 Her system than in HR 4049. This is the case, in particular, for a homogeneous dust composition throughout the disk, as is assumed here.

8.5. The scattered light

Finally, we discuss the main result of this paper, which is the inability of our disk models to reproduce the detected optical scattered light. Despite the extensive grid of calculated models, and the two dust compositions that were explored, not a single model produces sufficient scattered light at optical wavelengths. Moreover, those models that approach the required optical and resolved flux fraction do so at the cost of a too large scattering surface area and/or too high near-IR fluxes.

Considering the single-scattering albedos in Fig. 8, it becomes clear that another change in dust composition will not increase the scattered light by much. Doubling the iron abundance to 30% in mass fraction (an extreme situation) only slightly affects the albedo. One might expect that removing the smallest grains (≤0.03 μm) from the population, and including more grains in the ~0.03–0.3 μm size range by adapting the grain size power-law index, might increase the amount of scattered light. In contrast to the former, the latter grains have a high albedo in the optical wavelength range, and do not yet have a scattering phase function that is strongly forward-peaked. Mulders et al. (2013) showed that this strong forward peak in the scattering phase function of grains ≥0.3 μm is the cause for the faintness of protoplanetary disks in optical scattered light. We therefore recomputed our best-fit model, now excluding the smallest grains, and found no noticeable difference in the temperature structure nor in the predicted observables of our model. Although the smallest grains have high extinction coefficients, their density (see Fig. 10) is very low and they thus do not dominate the opacity in this model.

The small deficiencies of our best-fit Fe model with respect to the SED, mentioned in Sect. 5.1, point to a lack of small grains, and hence a power-law index that has too shallow a grain size distribution. Increasing the number of small grains is, however, prohibited by the constraints set in the H band. Our best amorphous carbon model partly resolves some of these issues, but still overestimates the H band flux and underpredicts the far-IR flux. Moreover, the origin of the carbon would be difficult to explain in this system. Finally, the visibilities are not perfectly fitted by either model, as an even flatter slope seems to be required at the intermediate-to-long spatial frequencies.

We conclude that no reasonable disk model can match all of the observations, because two distinct geometric structures are needed to accommodate the near-IR and optical constraints simultaneously. The structure of our best models is too flat at the inner rim to scatter a sufficient amount of the optical stellar flux at this radius into our line of sight, which is needed to fit our optical observations. In contrast, this flattened inner rim is needed to fit the near-IR data, and in particular the visibilities. So to explain the deficiencies, we propose that the basic geometry needs to be adapted: an outflow in optical scattered light should be included in the model of the system in addition to the circumbinary disk. As discussed in Paper I, there is some indirect evidence from polarimetric data (Joshi et al. 1987) that 89 Her may also have a low-density component of small dust grains that can be associated with a bipolar structure.

8.6. The global picture

Several questions remain about our interpretation of this system. The CO hourglass outflow is detected on a very different scale (~10″, Bujarrabal et al. 2007) than the scales traced by our interferometric observations. Moreover, the same authors estimated the lifetime of this hourglass at ~3500 yr. During this time the dynamics of the inner system may have changed. The Hα line shows a clear P Cygni type of profile (Waters et al. 1993), but the properties of the current outflow may be very different (e.g., more jetlike) from the CO hourglass. Based on our current modeling of the existing observations, we cannot conclude whether the optically scattering dust is dynamically linked to the hourglass, nor to a speculative more recent jetlike outflow.

The driving mechanism of the CO hourglass is unknown, so there are no independent constraints on the density structure nor on the grain properties of a potential dust component in it. Direct dust condensation in a wind seems unlikely given the effective temperature (6550 K) of the star, but dust grains from the circumbinary disk may be trapped in an accretion flow that is subsequently blown away in the polar directions. The launching mechanism of the gas and dust (i.e., the exact path followed by either) remains an open question in this framework (see, e.g., the model calculations of Bans & Königl 2012). The smoothness of the inner disk rim as found in this work, may be an indication for such a scenario where ongoing accretion is linked to the formation of the outflow. More post-AGB disks should be studied with near-IR interferometry to investigate a possible link between the inner disk shape and the presence of an outflow in the system.

Additional quantitative models that include the disk as well as a dust component in the (hourglass) outflow are needed in combination with more constraining interferometric observations. The parameter values of our best-fit disk models may be affected by the presence of dust in the outflow. It cannot be excluded that the outflow also contributes near- and/or mid-IR flux, which was neglected in the current framework. It might even be expected that this additional structure of low optical depth in/above the system presents a natural way of smoothing the intensity profile when viewed face-on. With this additional structural component, however, many new free parameters would be introduced in the model, in addition to the disk parameter space that needs to be reevaluated simultaneously. Therefore, degeneracies can be expected that the current data are incapable of resolving. Given the correspondence between the SED and the albedo of 0.1 μm grains (see Fig. 8), we nevertheless predict that an outflow consisting of such grains, and covering a large volume but with a low optical depth, might give rise to a significant flux of optical scattered light without producing too much IR thermal emission. It is, however, beyond the scope of this paper to explore such models.

9. Conclusions

In this paper we have computed an extensive grid of radiative transfer models of passive disks in hydrostatic equilibrium, using a grain size distribution that includes large grains and a consistent treatment of grain-size-dependent dust settling. Such models very successfully match the energetics and spatial distribution of the circumstellar material in the 89 Herculis system, as probed by the full IR SED and our IR interferometric observables. We have extended the SED with precise SPIRE flux measurements in the 200–600 μm wavelength range, thereby putting strong constraints on the size distribution of the dust grains. Our results confirm that the binary in 89 Herculis is surrounded by a compact and stable circumbinary disk in Keplerian rotation, wherein large dust grains form and settle to the midplane as in our best models. We also showed that no firm conclusions can be drawn on the composition of the dust, based on the observables used in this paper, without taking the full disk structure into account.

In the optical, however, only about half of the required circumstellar scattered flux is predicted by our disk models, albeit at the correct angular size scale. Although the parameter space is vast, and not all avenues in terms of disk models have been investigated in full detail (see, e.g., Turner et al. 2014, where the possible role of magnetic fields is investigated in the case of disks around Herbig stars), our results show that to resolve the remaining discrepancies, an additional component with some scattering material in our direct line-of-sight to the central binary is probably needed in the 89 Her system. The very smooth emission profile at near-IR and the compactness at optical wavelengths are difficult to reconcile within a geometrical framework that only consists of a puffed-up circumbinary disk in hydrostatic equilibrium.

Whether the dusty material is in the (detected molecular) bipolar outflow, in a halo near the inner disk rim, or in yet another structure, is unclear at the moment. Krijt & Dominik (2011) suggested that halos could form near the inner edge of protoplanetary disks by dynamical interactions and collisions between planetesimals. In light of our speculative detection of a “gap” in the MIDI data of 89 Her, such a scenario cannot be excluded a priori. It seems more likely, however, that the dust is present in the bipolar outflow, given that dust is detected in a variety of outflows around other post-AGB objects. The central question that remains relates to the origin of the outflow and its dust component. The grains in the outflow need to be ≥0.1 μm to efficiently scatter up to wavelengths of ~1 μm. The details of dust condensation in the close neighborhood of oxygen-rich (post-)AGB stars are not well understood, so the origin of these grains remains to be demonstrated. The mass reservoir from which the gas component in the outflow originates is not firmly identified either. The mass could either be lost by the central star(s) or by the inner disk. More (precise) observations are required to make an exploration of these complex model geometries fruitful.

We stress that our best-fit models of the inner geometry as listed in Table 1 were determined under the assumption that there is only a binary with a circumbinary disk in the system! On the other hand, we also underline that even though the outflow likely contributes scattered light at optical wavelengths, the difference in size distribution between the grains in the disk and in the outflow means that the (far-)IR emission is definitely dominated by the disk.

The complex inner geometry, the large scattering component, and the unknown physical connection between the CO hourglass outflow at large scales and the inner disk, indicates that multi-wavelength interferometric imaging is likely needed to constrain a global model of the 89 Her system. More clues on the 3D position of the scattering material can also come from stronger observational constraints on the wavelength dependence of the optical scattering, and from spatially resolving the system in different polarization states, since the various geometries can be expected to result in different polarimetric signatures.

Acknowledgments

This research has made use of the SIMBAD database, operated at the CDS, Strasbourg, France. This research has made use of NASA’s Astrophysics Data System. The Palomar Testbed Interferometer was operated by the NASA Exoplanet Science Institute and the PTI collaboration. It was developed by the Jet Propoulsion Laboratory, California Institute of Technology with funding provided from the National Aeronautics and Space Administration. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research. We also thank B. Acke for his input and the useful discussions about this work. Finally, we are grateful for the constructive comments and questions provided by the referee.

References

All Tables

All Figures

Fig. 1 Three baselines of MIDI data overplotted (at a single continuum wavelength, indicated in cyan) with all double power-law MCMax models computed in Sect. 5. The dotted blue line, “Best 1” in Table 1, is the best single power-law MCMax model in terms of fitting the near-IR visibilities (see Figs. 4 and 5). The dashed red line corresponds to the best-fitting double power-law model with the metallic Fe opacity source (“Best 2” in Table 1). The full green line shows the best-fitting double power-law model when metallic Fe is replaced by amorphous C (“Best alt.” in Table 1). See Sects. 4.2, 5.1, and 6 for a discussion of the models. The inset shows the data at the 30 m baselines as a function of wavelength, with the same double power-law models overplotted as in the main panel.

In the text

	Fig. 2 PTI absolute fluxes (in the insets) and spectral shapes over the H (upper panel) and K (lower panel) bands. The data are shown in blue. The dashed red and full green lines correspond to the same models as in Fig. 1.
In the text

Fig. 3 Spectral energy distribution of 89 Her with the same three models as in Figs. 1 and 2 overplotted. The dot-dashed gray line shows the input Kurucz model, here reddened with an interstellar contribution of E(B − V) = 0.07 (see Paper I). The black and gray data points are the “stellar” V,R,H, and K band fluxes as determined in Paper I. Open symbols represent actual measured fluxes in the passbands that are listed in the legend. The full black line is the ISO spectrum. The dotted blue line, “Best 1” in Table 1, clearly does not agree with the observed SED. A better, but still not excellent, fit to the SED can be obtained with an inner radius twice as large, in combination with p = 1.5, but this results in first-lobe near-IR visibilities that are too small (see the text). The dashed red and full green line correspond to the best-fitting double power-law models with metallic Fe and amorphous C opacity sources, respectively (see Fig. 1). The MCMax models have been reddened similarly as the Kurucz model.

In the text

Fig. 4 H band interferometric data. The left and right panels show the H band visibility amplitudes and squared visibilities observed by the various instruments. The visibility amplitudes, as produced by the data reduction software packages of the CLIMB and CLASSIC instruments, are not squared in order to preserve the original error statistics. The data are the same as in Paper I. The same three models are shown as in Figs. 1 and 3 and as listed in Table 1.

In the text

	Fig. 5 Same as Fig. 4 but for the K band.
In the text

	Fig. 6 Mass absorption (left) and total extinction (right) coefficients of the mixtures used in this work. Our “standard” mixture (full lines) contains 15% metallic iron, while the alternative mixture (dashed lines) has 15% amorphous carbon instead. The latter composition is exactly the same as that used by Mulders & Dominik (2012). Only six grain sizes are shown for clarity.
In the text

	Fig. 7 NPOI and VEGA squared visibilities at 673 nm. Our complete grid of double-power-law MCMax models is also displayed (see Sect. 5). The best-fit models listed in Table 1 are shown with the same color coding as in Figs. 1–5.
In the text

	Fig. 8 Albedo of the smallest three particle sizes and for four different compositions: silicates + amorphous carbon (dashed lines), only silicates (dotted lines), silicates + 15% iron (full lines), and silicates + 30% iron (dot-dashed lines).
In the text

Fig. 9 Temperature structure of our best-fit models. The upper panels show maps of the temperature, with the radial and vertical τ = 1 surfaces at 0.55 μm overlaid as solid and dashed lines, respectively. The middle and lower panels display the midplane pressure scale height and integrated dust density, respectively, as a function of the disk’s radius. On the left and right, the best Fe and amorphous C containing models are shown, respectively.

In the text

Fig. 10 Density structure of our best-fit models. From top to bottom, maps of the gas density, the dust density, the dust density for particles of ~0.06 μm, and the dust density for particles of ~5 mm are shown. On all panels the radial and vertical τ = 1 surfaces at 0.55 μm are overlaid as well, with solid and dashed lines, respectively. The left and right panels contain the best Fe and amorphous C containing models, respectively.

In the text