GAS DISTRIBUTION, KINEMATICS, AND EXCITATION STRUCTURE IN THE DISKS AROUND THE CLASSICAL Be STARS β CANIS MINORIS AND ζ TAURI^*

In order to obtain a model for the photospheric emission of β CMi, we searched for a consistent set of stellar parameters, taking the evidence for near-critical rotation into account (Saio et al. 2007). For this purpose, we employed our rapid rotator code (Monnier et al. 2007; Che et al. 2011), which simulates the stellar oblateness and surface temperature distribution using the modified von Zeipel theorem (gravity darkening coefficient = 0.188), and computes model images in the continuum and in the photospheric Brγ absorption line.

In order to constrain the stellar parameters, we computed a small parameter grid, in which we systematically varied the inclination angle i and the fractional angular velocity ω/ω_crit (where ω_crit is the critical angular velocity). For a given inclination angle and fractional angular velocity, we start our iterative process by assuming a stellar mass M_⋆, from which we compute the polar radius R_pole using the rotation velocity vsin i = 244 ± 6 km s⁻¹ (Yudin 2001). The polar temperature T_pole is derived from the V-band magnitude (V = 2.89). Any disk continuum flux in the optical band will tend to dilute the photospheric spectrum causing absorption lines to appear weaker than expected. We checked for such line reduction in blue spectra of β CMi obtained by Grundstrom (2007). Figure 5 shows the observed profiles (solid line) of He i 447.1 nm and Mg ii 448.1 nm based upon the average of 11 spectra made between 2005 and 2009 with the Kitt Peak National Observatory Coude Feed Telescope (spectral resolution R = 12, 500). The observed spectrum is compared with a theoretical flux spectrum (dashed line) from the grid of LTE models by Rodríguez-Merino et al. (2005). We adopted average parameters for the visible hemisphere of the star of T_eff = 11, 800 K, log g = 3.8, and Vsin i = 230 km s⁻¹ (Frémat et al. 2005), and the model spectrum was convolved with a simple rotational broadening function for a linear limb darkening coefficient of = 0.42 (Wade & Rucinski 1985). Note that this approach assumes a spherical star, ignores gravity darkening, and neglects changes in the local intensity spectrum with orientation, but these simplifications are acceptable for our purpose of checking for systematic line depth differences. We see that the match of the Mg ii λ4481 profile is satisfactory, but the wings of the He i λ4471 line appear too high, probably due to low level emission in this transition. Adding a 20% disk flux contribution and renormalizing the spectrum yields a diluted version (dotted line), which appears to be much too weak compared with the observed spectrum. This comparison suggests that the disk continuum emission in the blue range is negligible, and thus we will ignore any disk contribution to the optical flux in the following discussion.

Zoom In Zoom Out Reset image size

With the derived polar radius and bolometric luminosity, we are able to locate the stellar position on the HR diagram, after correcting for the rotational effect (Che et al. 2011). The stellar mass from the HR diagram is used as initial value for the next iteration step, until convergence between the assumed mass and the mass estimated from the HR diagram is reached.

For each parameter combination, our model provides the rotation period, apparent luminosity L^app, and apparent effective temperature T^app_eff, as listed in Table 2. Comparing the model effective temperature and luminosity with the observational constraints for β CMi (T_eff = 12, 050 K; L = 195 ± 60 L_☉; Saio et al. 2007), we find good agreement for models A, B, and D. Due to its consistency with the inclination and effective temperature (Sections 4.2 and 4.5), we favor model D, suggesting i = 40°, ω/ω_crit = 0.99, T_pole = 14, 550, vsin i = 244 km s⁻¹, T^app_eff = 12 871 K, and a polar and equatorial radius of R_pole = 2.95 R_☉ = 0.26 mas and R_eq = 4.10 R_☉ = 0.36 mas, respectively. We plot the values corresponding to our model grid in the HR diagram shown in Figure 6 and compare it to the evolutionary tracks from Yi et al. (2003) and Demarque et al. (2004), yielding an evolutionary age of ∼0.1–0.15 Gyr.

Zoom In Zoom Out Reset image size

Obviously, an important input parameter for our modeling procedure is the rotation velocity vsin i, for which we use an average value from the literature. Townsend et al. (2004) argued that most vsin i measurements on Be stars might systematically underestimate the true projected rotation value due to the effect of gravity darkening. In order to test this scenario, we have artificially increased the measured vsin i-value by 10% and repeated our modeling procedure. We find that the increase in vsin i results in a significant increase in the stellar mass, apparent effective temperature, and apparent luminosity, making the model prediction much less consistent with the observed values. Therefore, we suggest that in the case of β CMi, the vsin i-measurements are not significantly biased by gravity darkening, which is likely a result of the intermediate inclination angle of β CMi (the effect would be stronger for an equator-on viewing angle) and the use of our moderate gravity darkening law (T_eff∝g^0.188_eff, Che et al. 2011), which is likely more realistic than the classical von Zeipel darkening coefficient (T_eff∝g^0.25_eff) adopted by Townsend et al. (2004). Clearly, the accurate determination of vsin i remains a fundamental problem for constraining the stellar parameters of near-critical rotating stars. Therefore, we are currently working on incorporating the computation of model spectra in our modeling procedure, which will then be fitted to observed spectra together with the photometric and interferometric constraints (Che et al., in preparation).

The modeling provides spectra as well as model images in the continuum emission and in hydrogen photospheric absorption lines, which we will use as a representation of the stellar brightness distribution for our modeling of the disk in the following sections.

With projected baseline lengths up to B = 314 m, our CHARA/MIRC interferometric observations allow us to constrain the disk geometry in eight spectral channels in the H-band with an effective resolution of λ/2B = 0.5 mas. A commonly applied modeling approach is to approximate the disk emission with a Gaussian intensity profile, which is superposed on a uniform disk representing the stellar surface (e.g., Tycner et al. 2005; Meilland et al. 2009; Pott et al. 2010; Schaefer et al. 2010). A caveat of this approach is that a significant part of the disk emission is concentrated at radii r < R_⋆, and therefore wrongly attributed to the disk instead of the stellar flux. Accordingly, these models will systematically overestimate the disk emission with respect to the stellar emission (F_disk/F_⋆), resulting in inconsistencies with spectral energy distribution fitting results. In order to avoid these biases, we employ an elliptical Gaussian model, where the radial intensity profile at r > R_⋆ is given by a half-Gaussian, resulting in a proper estimation of the F_disk/F_⋆ ratio. As alternative model, we considered an inclined ring model, where the radial intensity profile is given by a Gaussian centered at radius a from the star, with a fixed fractional width of 25% (see Monnier et al. 2006a for details). Based on radiative transfer simulations of the brightness distribution in classical Be star disks (e.g., Waters 1986), we consider the Gaussian model a better representation of the expected emission profile in a Be disk, while the ring model is more suited for comparison with results from the literature (e.g., Meilland et al. 2006).

For all models, the stellar photosphere is represented by the aforementioned rapid rotator model (Section 4.1) and the position angle θ of the stellar equator and the disk major axis are aligned. Besides the disk major axis a (ring major radius or Gaussian half-width half maximum), the PA θ, and the inclination i, we also treat the flux ratio between the stellar and disk component (F_disk/F_⋆)_H as a free parameter. The model was fitted to the MIRC squared visibilities, closure phases, and triple amplitudes using a least-square fitting procedure, resulting in the best-fit model shown in Figure 7. The best-fit model parameters are displayed in Table 3, including 1σ-errors, which we have determined using a boot strapping procedure. We also tested whether the agreement can be improved using a skewed ring model (Monnier et al. 2006a), but find that introducing any disk asymmetry only marginally improves the fit. Our model includes minor asymmetries due to the brightened polar region in our rapid rotator model, which results in small closure phases (≲ 12), consistent with the measurement (Figure 2).

Zoom In Zoom Out Reset image size

The determined disk position angles of 1392 (ring model) and 1393 (Gaussian model) are in excellent agreement with the gas disk rotation axis determined with our VLTI/AMBER photocenter analysis (1400 ± 17).

Using the detailed information about the H-band geometry obtained with our extensive MIRC data set, we then fitted our Gaussian model to the CLIMB data to determine the K-band continuum geometry. Given the lower amount of observational constraints for the K-band, we treated only the disk-to-star flux ratio (F_disk/F_⋆)_K as a free parameters and kept the remaining parameters fixed. The resulting best-fit visibility curves are shown in Figure 7. With (F_disk/F_⋆)_K = 0.25^+0.09_−0.08, we do not find any significant deviations between the K-band and H-band flux ratio.

With the currently achievable accuracy, DP measurements can already reveal photocenter displacements two to three orders smaller than the formal angular resolution (λ/2B). For instance, our 2009 December 31 AMBER observations exhibit a DP accuracy of ∼08 (standard deviation over all continuum spectral channels), corresponding to a photocenter displacement of ∼8 microarcsecond (0.008 mas) on the employed 115 m baseline. Given that this is much smaller than the equatorial stellar diameter (∼720 microarcsecond), it is important to investigate whether the measured DP signatures might also contain contributions from the photospheric Brγ absorption line, which is tracing the stellar rotation and is underlying the Brγ-line emission.

In order to simulate the effect of stellar rotation on our interferometric observables, we employ our rapid rotator code and compute the stellar surface brightness distribution around the Brγ-line with a similar resolution as our AMBER observations (Figure 8, top). For this, we align the stellar rotation axis with the measured rotation axis (θ = 2275, Section 3). Brγ absorption on one side of the photosphere will shift the photocenter toward the opposite direction. Taking this into account, we adjust the orientation of our model photosphere, so that the Brγ absorption in the blue-shifted line wing matches the measured direction of the photocenter displacement at red-shifted wavelengths. In addition to the photospheric emission, we include disk emission in our model, assuming (F_disk/F_⋆)_K = 0.25, as determined in Section 4.2.

Zoom In Zoom Out Reset image size

From the model, we compute the corresponding DPs and find signatures ϕ ≲ 05 on all baselines (Figure 8, bottom). This result reflects the fact that the fraction of the photospheric absorption to the total K-band emission is rather small (<8%) and that the absorption appears as a rather diffuse structure on the extended stellar surface. We conclude that the photocenter signatures due to stellar rotation (in photospheric absorption lines) are a minor effect compared to the rotation signatures of the disk (in line emission), but nevertheless result in DP signatures comparable to the currently achievable DP accuracy and are therefore important, in particular for more detailed future studies.

Position–velocity (p–v) diagrams provide a powerful tool for the interpretation of a disk velocity field and are commonly employed in radio astronomy to derive the rotation profile of circumnuclear galactic disks or protostellar disks using, for instance, maser (e.g., Miyoshi et al. 1995; Pestalozzi et al. 2009) or molecular emission tracers (e.g., Sofue & Rubin 2001; Isella et al. 2007). Using our photocenter displacement measurements, we can construct an equivalent diagram from our VLTI interferometric data by measuring the length of the continuum-corrected photocenter displacement vectors projected on the disk plane θ = 1400. Performing this projection on the disk plane also allows us to avoid a potential bias due to opacity effects, since these effects would move the photocenter only perpendicular to the disk plane (Section 3). The resulting p–v diagram (Figure 9) shows a symmetric rotation curve, which we interpret using the simple model of a thin Keplerian-rotating disk (Weintroub et al. 2008). The disk extends from an inner (R_in) to an outer (R_out) radius and the line-of-sight (LOS) velocity v of the emission element located at radius r and at angle ϑ in the disk plane is given by

$$\begin{equation} v_{\rm kep}(r,\vartheta) = \sqrt{\frac{GM_{\star }}{r}} \sin \vartheta \cdot \cos i, \end{equation} \tag{ 4 }$$

where ϑ is measured against the LOS and G is the gravitational constant. Considering only the emission from the outermost disk annulus (i.e., at R_out and ϑ = −π... + π) will result in a straight line in the p–v diagram (line A-B in Figure 9). Adding the emission from the remaining disk annuli (R_in to R_out) will result in a characteristic "bowtie"-shaped filled region in the p–v diagram, such as commonly observed in CO imaging observations (Sofue & Rubin 2001; Isella et al. 2007). However, astrometric observations (such as radio masers or photocenter displacements) trace the light barycenter of each annulus, corresponding to the B-C and A-D curves in Figure 9 (for more details see Weintroub et al. 2008). We fit this simple model to the β CMi p–v diagram assuming a fixed inclination of i = 40° (as determined in Section 4.2), which yields best agreement for M_⋆ = 3.6 ± 0.3 M_☉, R_out = 2.4 ± 0.2 mas, and R_in < 1.0 mas (solid curve in Figure 9).

Zoom In Zoom Out Reset image size

By reducing the information content to purely astrometric data, the p–v diagram analysis method presented in the last section provides a very intuitive and powerful method to constrain the gas velocity field in the observed line tracer. In this section, we will use a more sophisticated model in order to make full use of the rich information contained in our spectro-interferometric observations, including spectra, wavelength-differential visibilities, DPs, and CPs. The Brγ-emitting gas is assumed to be located in a thin disk plane, which is justified due to the expected small opening angle of Be star disks (e.g., Bjorkman & Cassinelli 1993) and the intermediate inclination angle under which β CMi is observed (Section 4.2).

The gas kinematics is parameterized with a rotation profile $|\vec{v}(r)| = f_{\rm kep}(R_{\rm ref}) \cdot (r/R_{\rm ref})^{\beta }$, where β = −1 for rotation with constant angular momentum, β = −0.5 for Keplerian rotation, β = 0 for constant rotation, and β = +1 for solid body rotation (e.g., Stee 1996). $f_{\rm kep}=|\vec{v}(R_{\rm ref})| / |\vec{v}_{\rm kep}(R_{\rm ref})|$ is the orbital velocity at the reference radius R_ref = 1 AU, expressed as fraction of the Keplerian velocity $|\vec{v}_{{\rm kep}}(r)| = (G M_{\star }/r)^{-1/2}$. In the model, the line emission extends from the stellar surface to an outer radius R_out with a radial power-law intensity profile I_l(r)∝r^q. To include thermal line broadening in our kinematic model, we adopt a constant radial gas temperature of 0.6 · T_eff (=7860 K for β CMi), as suggested by the radiative transfer modeling from Carciofi & Bjorkman (2006).

We include both line and continuum emission in our model and compute the interferometric observables for our given VLTI array configurations and the covered wavelength channels. The wavelength-dependent model visibilities and phases are then fitted to the VLTI interferometric data using a reduced χ²_r goodness-of-fit estimator (see Kraus et al. 2009 for a definition), including our flux, visibility, DP, and CP constraints.

For the continuum emission, we assume the geometry determined with our CHARA observations (Section 4.2) and use this model to renormalize the AMBER continuum visibilities. In order to incorporate the underlying photospheric Brγ absorption, we include the photosphere model discussed in Section 4.3, although, as discussed above, the influence on the differential phase (≲ 05 on all baselines) is still within the measurement uncertainties.

The disk position angle is fixed to θ = 1400, as determined by our model-independent photocenter analysis (Section 3). The remaining six parameters in our modeling are the outer disk radius R_out, the inclination i, the stellar mass M_⋆, the radial intensity power-law index q, the disk rotation index β, and the velocity at the reference radius R_ref (expressed as fraction of the Keplerian velocity, f_kep). In a first step, we fix the velocity profile to Keplerian rotation (β = −0.5, f_kep = 1) and vary the remaining four parameters systematically on a grid, yielding the best-fit model shown in Figure 10. The corresponding parameters and uncertainties are listed in Column 4 of Table 3. In a second step, we test also non-Keplerian velocity fields, yielding the best-fit values listed in Column 5. The resulting χ²-surfaces are shown in Figures 11 and 12 and the uncertainties have been derived using the bootstrapping technique.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our model fits show that the intriguing phase inversion observed at the line center on our longest VLTI baselines (Figure 10, 3rd row) can be explained with the phase jumps in the Fourier phase crossing a visibility null. These phase jumps appear in the same spectral channels where we measure visibility minima in the W-shaped visibility profile, indicating that the visibility function of the pure line-emitting geometry transits here from the first to the second visibility lobe (Figure 10, 2nd and 3rd row). At the same time, the continuum emission is only marginally resolved, which results in the measured composite line+continuum visibility with a rather high contrast of ∼0.7. The phase jumps and visibility minimums are a basic property of the Fourier transform of strongly resolved objects and are reproduced very naturally and without fine-tuning from our kinematical modeling. Therefore, our results do not support the idea outlined by Štefl et al. (2011) that these features might indicate secondary dynamical effects or polar mass outflows.

As best-fit value for the radial intensity index, we yield q = −1.6 ± 0.2. Isothermal viscous decretion disk models, such as discussed in Bjorkman & Carciofi (2005), predict a radial disk surface density law Σ(r)∝r⁻², corresponding to q = −2 in the optically thin case. More realistic non-LTE disk models suggest a more shallow surface density profile (Σ(r)∝r^{−2... − 1}) due to a steep temperature drop in the inner few stellar radii of the disk (Carciofi & Bjorkman 2008). Therefore, we conclude that our measured intensity profile is well consistent with these models.

The most significant deviation of our simple kinematical model from the data is in the precise shape of the Brγ-spectrum. For instance, the measured spectra show a weak emission component in the line center (Figure 10), which is not reproduced by the model and which might be related either to the presence of a uniformly distributed low-velocity gas halo or, more likely, to radiative transfer effects.

Another deviation between the model and the observation concerns the weak V/R-asymmetry which can be observed in our VLTI/AMBER Brγ spectra at both epochs (Figure 3, top row) and which is not reproduced by our axialsymmetric kinematical model. Within the 113 days covered by our AMBER observations, no changes in the V/R asymmetry could be observed, which is consistent with the conclusion by Tycner et al. (2005) and Jones et al. (2011) that the Hα profile does not show significant long-term variability. It is interesting to compare the average photocenter displacement between the blue- and red-shifted line emission (2.5 mas) with the characteristic size of the Hα-emitting region (2.13 mas, Tycner et al. 2005), which suggests that Brγ emerges from a similar spatial region as Hα and a much more extended region than the H- and K-band continuum emission (R_Brγ ∼ R_Hα > R_cont, Figure 3), in agreement with the prediction from Carciofi (2010). However, a detailed comparison is difficult due to the non-simultaneity of the different observations.

Both in Brγ and in the Pf14-Pf22 lines, we detect a double-peaked line profile and clear rotation signatures in the differential phases. All line transitions exhibit a photocenter displacement with a stronger amplitude in the southeastern (blue-shifted) than in the northwestern (red-shifted) lobe (Figure 4, middle panel). Such an asymmetric displacement is consistent with the presence of a one-armed oscillation in the disk (Štefl et al. 2009; Carciofi et al. 2009). CHARA/MIRC observations by Schaefer et al. (2010) showed that the one-armed oscillation pattern can also be observed as an asymmetry in the H-band continuum emission. Using multi-epoch data, they find that the position angle of the asymmetry is correlated with the spectroscopic variability and precesses around the star with the Hα V/R period. Our observation (2010 January 1) adds an additional epoch for this analysis at the V/R phase 0.484 (assuming maximum phase at JD = 2454505.0 and a period of 1429 days, Schaefer et al. 2010). We find that the spiral density maximum is located toward the southeast of the central star and the measured PA is consistent with the sinusoidal PA modulation proposed by Schaefer et al. (2010), as shown in Figure 13. Based on our independent calibration of the DP-sign (see Section 3), we find that the approaching (blue-shifted) part of the ζ Tau disk is located southeast of the star, which corresponds to an opposite rotation sense than the one adopted by Štefl et al. (2009).

Zoom In Zoom Out Reset image size

The reported polarization angle for ζ Tau (329–331, Quirrenbach et al. 1997; 35° ± 4°, Ghosh et al. 1999; 32°–335, McDavid 1999) is in excellent agreement with the rotation axis position angle (35.8 ± 2.0) found by our spectro-interferometric observations.

In order to test whether the measured differences in the stellocentric emitting radius of the Brγ and Pfund transitions are consistent with the expected excitation structure in the disk, we construct a simple radiative transfer model assuming local thermodynamic equilibrium (LTE). We assume that the line-emitting material is located in an equatorial disk which extends outward from the stellar radius with a constant vertical density per unit volume and a half-opening angle of Θ = 5°. Then, we integrate for each radius r the optical depth τ_ν in vertical direction, assuming hydrogen under LTE conditions and an isothermal temperature distributions. The number density for the different excitation levels and ionization stages is computed using the Saha and Boltzmann equation assuming a Gaussian line profile with thermal line broadening (Wilson et al. 2009). The Einstein coefficients for the different hydrogen transitions are estimated using the series expansion published by Omidvar & McAllister (1995). The radial temperature and surface density profile are parameterized with T(r) = T_eff(r/R_⋆)^−0.5 (Stee & de Araujo 1994) and Σ(r) = Σ₀(r/R_⋆)⁻² with Σ₀ = 2.1 g cm⁻² (Carciofi et al. 2009). For the stellar temperature, equatorial stellar radius, and distance, we assume T_eff = 19 370 K, R_⋆ = 7.7 R_☉, and d = 126 pc, respectively (Carciofi et al. 2009).

Using the radiative transfer equation for LTE conditions ($I_{\nu }(r)=B_{\nu }(T) (1-e^{-\tau _{\nu }(r)})$), we compute the emitted intensity per unit area as function of radius r in the disk, where B_ν is the Planck spectrum for temperature T at frequency ν. From the radial intensity profiles (Figure 14) we compute the centroid of the brightness distribution, which is then compared to the stellocentric emission radius measured in the different line transitions. Given that our model assumes a simplified vertical density structure and does not include inclination effects, we do not aim to match the absolute sizes of the emitting region in all line transitions, but focus instead on the relative sizes. For this step, we normalize both the measured and the model photocenter offsets relative to Brγ. The comparison between the Hα, Brγ, and Pf14-17 relative sizes and our LTE model is shown in Figure 15. We find that we can reproduce important observational features, in particular that

• 1.  
  Brγ originates from a similar spatial region in the disk as Hα (R_Brγ ∼ R_Hα). This result is also in agreement with the predictions from more sophisticated non-LTE radiative transfer computations, such as made by Carciofi (2010). In order to better characterize the differences between the Brγ and Hα-emitting region, contemporaneous observations in these two wavelength bands will be required.
• 2.  
  Brγ originates from a more extended region than the Pfund lines (R_Brγ > R_Pf). Computing the weighted average of the measurements in the individual transitions, we yield $R_{\rm Pf}/R_{\rm {\rm Br}\gamma }=0.48 \pm 0.12$. For the Pfund lines, we are not aware of earlier predictions obtained with non-LTE radiative transfer codes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our observational results confirm the finding from Pott et al. (2010), which found R_Brγ > R_Pf on the classical Be star 48 Lib, and consolidates their suggestion that the measured size differences can already be explained with the expected optical depth differences between these line transitions. We encourage theoreticians and modelers to employ their more sophisticated non-LTE radiative transfer codes in order to test the influence of inclination and non-LTE effects, although, based on the results from Iwamatsu & Hirata (2008, e.g., their Figure 3), we expect no significant departure from LTE for the inner disk regions and the high transitions traced by our observations. Future multi-transition spectro-interferometric observations with improved uv-coverage might also measure the precise radial intensity profile in a model-independent fashion. Together with sophisticated radiative transfer simulations, these observations will reveal the excitation structure of the disk and constrain parameters such as the temperature profile and the vertical disk structure, which are currently difficult to access.