© ESO, 2015

1. Introduction

The gaseous circumstellar disks around classical Be stars are special universal laboratories. Although Be stars have been known for 150 years (Secchi 1866) and the idea of a flattened envelope has existed since 1931 (Struve 1931), consensus regarding the physical process governing the disk structure, the viscous decretion disk (VDD) model, has emerged only recently.

The VDD model, first introduced by Lee et al. (1991) and further developed by, e.g., Bjorkman (1997), Porter (1999), Okazaki (2001), and Bjorkman & Carciofi (2005), is now widely accepted as the best physical model for describing the circumstellar disks of Be stars. Among the growing evidence supporting the VDD model is the confirmation that the disks rotate in a Keplerian way (Meilland et al. 2007; Kraus et al. 2012; Wheelwright et al. 2012), allowing for the identification of viscosity as the mechanism that makes the disk grow: Viscous torques transfer angular momentum from the base of the disk outward, thus allowing the gas particles to reach progressively wider orbits. Decretion disks share most of the physical characteristics with accretion disks known, e.g., from protostars, with their name referring to the outward direction of the mass flow. The mass gets injected into the decretion disk by an as of yet not unambiguously defined mechanism, most probably a combination of nonradial pulsations, fast rotation, and possibly small-scale magnetic fields (Rivinius et al. 2013a).

Stringent tests of the VDD model focusing on reproducing the observables of individual targets are however still very few in number. The VDD model implemented in the Monte Carlo radiative transfer code HDUST (Carciofi & Bjorkman 2006, 2008) was used to reproduce the violet-to-red peak ratio (V/R) spectroscopic variability along with photometry, polarimetry, and near-infrared (IR) interferometry of ζ Tau, providing firm evidence that the V/R oscillations are an effect of one-armed density waves in the disk (Carciofi et al. 2009). Other successful HDUST applications include the visual light-curve modeling of the 2008 outburst of 28 CMa along with the first determination of the viscosity parameter (Carciofi et al. 2012), and the modeling of visual and IR spectral energy distribution (SED) and visual spectropolarimetry of δ Sco (Carciofi et al. 2006). Further studies of individual stars include those of Jones et al. (2008) and Tycner et al. (2008) in which the BEDISK code (Sigut & Jones 2007) was used. Circumstantial evidence supporting the VDD model is also provided, e.g., by Haubois et al. (2012), who were able to reproduce observed light curves of Be stars with dynamical models; by Silaj et al. (2010), who reproduced generally well the observed Hα profiles of 56 Be stars; and by Touhami et al. (2011), who were able to reproduce the statistical properties of the color excesses of a sample of 130 stars. The basic characteristics of the VDD model and HDUST code are described in Sect. 3. For an extensive review on classical Be stars and their disks, see Rivinius et al. (2013a).

When studying a complex astrophysical system, such as a star irradiating a surrounding circumstellar disk, it is crucial to combine different observational techniques, as each probes different aspects of the disk physics. The combination of radiative transfer and disk dynamics governs most observables, and their disentangling is the central goal of the modeling process. The most characteristic observables of Be stars are the often double-peaked hydrogen emission lines and the disk excess IR emission due to free-free and bound-free emission in the ionized circumstellar disk. Observed optical and near-IR radiation is partially linearly polarized as a result of scattering off free electrons in the inner disk. Medium- and high-resolution differential spectrointerferometry across an emission line, such as Brγ, also shows patterns in visibilities and phases that are characteristic of circumstellar disks rotating in a Keplerian way. The combination of spectroscopic, photometric, polarimetric, and interferometric observations is therefore particularly strong for the case of circumstellar disks, removing some of the degeneracies arising from fitting single observations.

The main goal of this study is to apply the VDD model to the case of the well-known Be star β CMi and to test how well it handles arguably the largest set of multiwavelength and multitechnique observations of a Be star so far, thereby broadening the sample of detailed VDD tests applied on individual targets. The data set includes a compiled SED covering the ultraviolet (UV) to radio wavelengths, high-resolution spectroscopy of the emission lines, polarimetry, and high-resolution spectrointerferometry (see Sect. 2).

Besides being nearby, bright and therefore often observed, β CMi was chosen as the target for this study for two additional reasons. Firstly, the disk of β CMi has been been present and stable for decades. Slettebak (1982) reported “no appreciable change” in line emission in the years 1950 to 1982 and the same holds for measurements from 1989 and 1991, which were published in Hanuschik et al. (1996). This implies that the disk does not go through major mass ejection and dissipation episodes often observed in many Be stars. The absence of high-amplitude V/R variations is also evidence of no large-scale density fluctuations throughout the disk. The uncertainties arising from combining observations from different epochs should therefore be minimized. Secondly, β CMi represents one of the first detailed VDD applications to a late-type Be star (B8Ve, Abt et al. 2002), thus enabling testing of the universality of the VDD model throughout emission line stars of the spectral class B. Moreover, it is not entirely clear whether late-type Be stars differ qualitatively from the early-type stars with respect to, e.g., the observed variability (which seems to be lower for later types) and the rotation rate (higher for later types, see Sects. 3.1 and 3.2 of Rivinius et al. 2013a).

Probing the outer reaches of the disk by analyzing radio observations is another main goal of our research. For cases nonedge-on, the SED fluxes can be interpreted as the superposition of the photospheric and disk contributions. The free-free mechanism dominates the disk opacity in the near-IR and longward, and increases with wavelength approximately as ∝λ². Consequently, in a typical Be disk the central parts are optically thick along the line of sight and act like a pseudophotosphere, which longward of a density-dependent λ₀ grows in size according to the approximate relation $\overline{\mathit{R}}\mathrm{\propto }{{}^{\mathrm{(}}\mathit{\lambda }\mathit{/}{{\mathit{\lambda }}_{\mathrm{0}}}^{\mathrm{)}}}^{\mathrm{0.41}}$, where $\overline{\mathit{R}}$ is the pseudophotosphere radius (Vieira et al. 2015). Therefore, different wavelengths probe distinct regions of the disk (see, for instance, Fig. 2 of Rivinius et al. 2013a). In this context, radio observations are of particular interest in studying the outer parts of the disk, and may eventually impose important constraints to the disk physical size. Radio observations were used to study the structure of the outer Be disks in the pioneering works of Taylor et al. (1990), Waters et al. (1991), and Dougherty et al. (1991b), however, since then no more radio studies of Be stars were published until very recently (see Sect. 2.2).

The binary nature of β CMi is still an open question. It was suggested as a spectroscopic binary by Jarad et al. (1989); however, their data cover less than two cycles (considering their period of 218.5 days) and the spectral lines in the interval 3700–4700 Å, which they used to determine the radial velocities, could be contaminated by disk emission. Their results for the radial velocities are therefore questionable and have not been confirmed. Recently, Folsom et al. (2015) reported weak V/R oscillations (V/R = 0.98–1.04) of the Hα emission line with the period of 182.83 days discovered by a spectroscopic monitoring campaign in the years 2000–2014. These variations may be interpreted as arising from disk asymmetries triggered by the tidal interaction between the disk and an unseen binary companion. In this case, the V/R oscillations would be phase-locked with the binary orbit and therefore the V/R period would be equal to the binary orbital period.

Binary interaction is the only known mechanism that can truncate the disk interior to the photoevaporation radius. The so-called truncation radius represents a position in the disk where the tidal torque caused by the companion is in balance with the viscous torque. At this radius the disk is not simply cut off, but rather separated into two density regimes: Within the truncation radius the radial density decrease dependence can become shallower than an equivalent, isolated disk, and outside the truncation radius it becomes steeper, while the magnitude of this effect depends on the viscosity parameter (Okazaki et al. 2002; Panoglou et al. 2015) and on the binary mass ratio (Sect. 5). If the reported V/R variations are indeed the effect of a binary orbit, radio observations might in principle be used to detect and quantify the truncation of the disk by the companion. In view of the possible change of the density structure with respect to steady viscous decretion caused by a binary truncation or a nonsteady decretion (Haubois et al. 2012), we investigate whether steady viscous decretion provides an adequate description or whether a more complicated density structure is necessary.

The last goal is to look for observable signatures of the disk nonisothermality, which is expected from theory. Be disks become nonisothermal because the inner, dense parts of the disk are optically thick to the photoionizing stellar radiation, causing an initial temperature decrease in the radial direction. As the density gets lower further from the star, the temperature rises to about 60% of the stellar effective temperature and the disk becomes roughly isothermal. The temperature structure is expected to have an effect on the density structure of the disk, which then affects the observables (Carciofi & Bjorkman 2008).

In Sect. 2 the observations are described, while Sect. 3 briefly describes the adopted modeling procedure. In Sect. 4 the model predictions are presented and with their help the observations are interpreted. The conclusions follow.

2. Observations

2.1. Optical and IR photometry

Standard Johnson UBVRI photometry was taken from the catalog of Ducati (2002). The UBVRI photometry was dereddened with the reddening curve of Cardelli et al. (1989) using the extinction coefficient R_V = 3.1 and the interstellar color excess E [B−V] = 0.01 adopted from the analysis of Dougherty et al. (1994). The overall reddening effects are found to be very weak (and negligible longward of 7000 Å), as is also indicated by the very low extinction toward β CMi (A_V = 0.05 according to Wheelwright et al. 2012).

IR photometry in filters JHKLM was adopted from Dougherty et al. (1991a). We also make use of the color-corrected IRAS measurements of β CMi published in Coté & Waters (1987). The remaining IR photometry was taken from the published measurements from the AKARI/IRC mid-IR all-sky Survey (Ishihara et al. 2010), from the Spitzer Space Telescope (SST; Su et al. 2006) and from the Wide-field Infrared Survey Explorer (WISE; Cutri et al. 2014). Finally, we use the N-band flux calibrated spectrum (R = 30) measured by the VLTI/MIDI interferometer. The compiled IR data set is listed in Table 1.

2.2. Radio measurements

Radio observations of Be stars have been very scarce throughout decades past. However, recent years saw a great technological improvement in the instruments allowing for precise measurements even of very low fluxes. We briefly discuss the historic measurements and then describe our own from recent years.

Already in the late eighties, several Be stars including β CMi were detected at 2 cm wavelength by the Very Large Array (VLA; Taylor et al. 1990) and at mm wavelengths by the James Clerk Maxwell Telescope (JCMT; Waters et al. 1991). We use the published JCMT and VLA data in our analysis (see the respective references for the description of the data). To our knowledge, these had been the only radio observations of Be stars up until a few years ago, when 28 CMa was observed by the Atacama Pathfinder Experiment (APEX) millimeter telescope (Štefl et al. 2011) and δ Sco was observed by APEX and by the Combined Array for Research in Millimeter-wave Astronomy (CARMA, Štefl et al. 2012). More recently, we used the APEX bolometer camera LABOCA (Siringo et al. 2009) and the CARMA array (Bock et al. 2006) in its photometric mode to observe β CMi. The sub-mm and radio observations are included in Table 1.

The APEX/LABOCA camera operates at the central frequency of 345 GHz, corresponding to the wavelength of 870 μm with a bandwidth of 150 μm. It consists of a 295-channel bolometer array laid out in concentric hexagons around a central channel. The sub-mm radiation is absorbed by a heat-sensitive semi-conductor, which measures temperature variations of a thin titanium film. Our measurement consists of five scans with a total integration time of 53.8 min. The data were reduced with the help of the free reduction software CRUSH¹ version 2.20-3 (Kovács 2008).

The CARMA observations were executed with the 15-antenna subarray (nine 6.1-m and six 10.4-m diameter antennas) in the D configuration, in which the distances between the individual antennas are 11–150 m. The central frequency of the observations was ~100 GHz (~3 mm) with the bandwidth of ~500 MHz. The total integration time including the calibrators was 2.60 hours. The data were reduced using the MIRIAD software package².

2.3. Spectroscopy

The spectral coverage of the compiled SED is extended into ultraviolet (UV) wavelengths by including spectra measured by the International Ultraviolet Explorer (IUE)³. We observed β CMi by IUE five times in high dispersion mode with the total spectral coverage of 1150–3350 Å. The spectra were averaged and then dereddened with the same approach as the optical photometry with the reddening curve updated for the near-UV according to O’Donnell (1994). Strong reddening is usually indicated by an interstellar absorption feature around 2200 Å. The fact that this feature is absent from the observed spectra is a further indication of low extinction toward β CMi and overall weak reddening effects.

In the optical region, the Johnson UBVRI photometry is complemented with several flux-calibrated low-resolution spectra covering the interval 3200–10 500 Å measured by the HPOL spectropolarimeter⁴ mounted on the Pine Bluff Observatory (PBO) 36′′ telescope. As the absolute flux calibration of the spectra is not reliable, however, only the shape of the spectra was used. The individual measurements were scaled to the V-band magnitude before they were averaged and dereddened.

The spectroscopic observations of hydrogen Balmer emission lines (R = 5000–30 000) in the 2001–2014 interval were carried out using the slit spectrograph mounted at the coudé focus of Ondřejov 2 m telescope in the Czech Republic and the Cassegrain (ECass) and MUSICOS spectrographs at the 1.6 m telescope at the Observatório Pico dos Dias (OPD) in Brazil. The gaps in the time coverage were filled in by several amateur spectra from the BeSS database⁵. The Balmer lines are complemented with Brγ profiles extracted from two sets of AMBER interferometric measurements (Sect. 2.5). Overall, the line profiles show only small-scale temporal variability, although some show a weak V/R assymetry, in line with the results of Folsom et al. (2015). In Fig. 1 we show the Ondřejov and AMBER line profiles, the averages of which were used for the modeling process.

Two spectra from the high-resolution (R = 48 000) fiber-fed spectrograph FEROS (Kaufer et al. 1999), mounted at the 2.2 m ESO telescope at the La Silla observatory, Chile and a set of spectroscopic measurements from the ESPaDOnS⁶ echelle spectrograph (R = 68 000), installed at the 3.6 m Canada-France-Hawaii Telescope (CFHT) in Hawaii, USA, allowed us to investigate the whole spectrum in the 3600–9200 Å range.

The information on the spectroscopic data set is summarized in Table 2.

2.4. Linear polarimetry

β CMi was observed by the HPOL Spectropolarimeter mounted on the PBO 36′′ telescope six times during the years 1991–2005. The first measurement was done with a Reticon dual array detector, while a CCD detector was used for the remaining scans. The wavelength coverage of the CCD observations is 3200–10 500 Å with a spectral bin size ~10 Å.

Recent broadband polarimetric measurements in the Johnson filters BVRI were obtained on five observing nights during the years 2011–2014 at 1.6 m telescope of the OPD observatory in Brazil (Table A.1). For the instrument description, observing setup and calibration and reduction procedures see Sect. 2 of Carciofi et al. (2007) and references therein.

As discussed in Sects. 2.1 and 2.3, the interstellar extinction toward β CMi is almost negligible, therefore, we perform no corrections for the interstellar polarization.

2.5. Interferometry

The interferometric observations were obtained at three interferometric facilities: Very Large Telescope Interferometer (VLTI) at the Paranal observatory in Chile, Center for High Angular Resolution Astronomy (CHARA) at the Mount Wilson Observatory in California, USA, and the Navy Precision Optical Interferometer (NPOI) located in Arizona, USA. The VLTI instruments include the 2-telescope beam combiners PRIMA FSU-A (Sahlmann et al. 2009) and MIDI (Leinert et al. 2003) and the 3-telescope beam combiner instrument AMBER (Petrov et al. 2007). We also observed β CMi for two nights (May 1–May 3, 2014) with the 4-telescope beam combiner PIONIER (Le Bouquin et al. 2011), but these measurements were mostly lost due to bad weather conditions. The CHARA instruments used are the 4 to 6-telescope beam combiner MIRC (Monnier et al. 2006) and a 3-beam combiner CLIMB (Sturmann et al. 2010). The AMBER, MIRC, and CLIMB data were already analyzed by Kraus et al. (2012) and part of the MIDI observations was published in Meilland et al. (2009). The observation logs for interferometric observations from VLTI and CHARA are listed in Table 3 and the full NPOI data set is given in Table A.2.

The former PRIMA (Phase-Referenced Imaging and Micro-arcsecond Astrometry) facility at VLTI consisted of two identical fringe sensor units (called FSU-A and FSU-B) but only FSU-A was publicly offered for observations. The FSU-A operates in the K-band between 2.0 μm and 2.5 μm. The light is chromatically dispersed over five spectral pixels. Three visibility measurements were obtained with pairs of 1.8 m Auxiliary Telescopes (ATs). For more details on the observing setup and the reduction procedure, see Sect. 4 of Müller et al. (2014).

The MIDI data consist of three measurements executed with pairs of 8.2 m Unit Telescopes (UTs) with spectral dispersion R = 30 across the N-band (8–13 μm). For more details on the data acquisition and reduction, see Meilland et al. (2009) and references therein.

The AMBER observations were obtained in the high spectral resolution mode (HR mode, R = 12 000) in the K-band spectral region centered around the hydrogen Brγ line (${\mathit{\lambda }}_{\mathrm{Br}\mathit{\gamma }}^{\mathrm{vacuum}}\mathrm{=}\mathrm{2.166078}\mathit{\mu }$m). The pointing from December 31, 2009 was obtained with three 8.2 m UTs, while for the observation from April 23, 2010 three 1.8 m ATs were used. For our purposes, we rereduced the data using the amdlib V3.0 data reduction software (Tatulli et al. 2007; Chelli et al. 2009). As absolute calibration of the AMBER data is not reliable, we normalized the continuum region and used the differential visibilities and phases for our analysis. As the measurement error for all data points we used the standard deviation computed from the continuum region. For more details on the observing setup and the reduction procedure see Sect. 2 of Kraus et al. (2012).

The observations from MIRC were obtained with low spectral dispersion, R = 35, in the H-band continuum. The light from four 1 m telescopes was combined with the baseline lengths reaching 330 m. The CLIMB observations include eight pointings in the K′-band continuum. The observing setup and the reduction process are detailed in Sect. 2 of Kraus et al. (2012). For the uv-plane coverage of the CHARA measurements as compared to AMBER coverage, see their Fig. 1.

The NPOI is a long-baseline, multielement, optical interferometer that can combine light beams from up to six elements simultaneously (Armstrong et al. 1998). Before the fringe contrast is recorded between the beams, the light combined at the beam combiner is dispersed using a prism onto a lenslet array connected by optical fibers to a series of photon-counting avalanche photodiodes, which results in data stream from 16 spectral channels covering the wavelength range of 560–870 nm. For B-type stars with an Hα emission line, it allows interference fringe contrast (in the form of squared visibility) to be recorded for a single 15-nm wide spectral channel that contains the entire emission line from the circumstellar region (we refer to this as Hα channel when the 656.28 nm line is centered in the spectral channel).

The NPOI observations of β CMi utilized in this study were obtained on 40 nights from November 2010 to April 2011. The entire data set resulted in 720 unique squared visibility measurements from the Hα channel. The processing of NPOI data was conducted using a custom suite of programs developed specifically for the instrument known as the Optical Interferometer Script Data Reduction package. The processing of the raw fringe data frames follows the procedures outlined in Hummel et al. (1998) with additional bias corrections using off-fringe measurements described in Hummel et al. (2003). Furthermore, in the case of Hα-emitting sources, where the signal of interest is confined to a single spectral channel per baseline, the additional step of removal of small channel-to-channel variations has also been conducted as described in Tycner et al. (2006) utilizing a calibrator star λ Gem (FK5 277, HR 2763). The last step of calibration of the Hα squared visibilities involved calibration with respect to the continuum channels, which are assumed to be represented by a uniform disk model with a known angular diameter (Tycner et al. 2003). In the case of β CMi, we adopted an angular diameter of the central star of 0.675 mas, which was obtained based on a Hipparcos distance of 49.6 pc (van Leeuwen 2007), the best-fit value for stellar polar radius (R_p) of 2.8 R_⊙ based on UV spectrum fitting (see Sect. 4.1.1) and a multiplicative factor of 1.29, which accounts for the average effect due to the rapid stellar rotation⁷.

3. Model description

We modeled β CMi as a combination of a fast-rotating central star surrounded by a VDD with the help of the computer code HDUST. The code HDUST uses the Monte Carlo method to solve the radiative transfer, radiative equilibrium, and statistical equilibrium in 3D geometry for arbitrary density and velocity distributions assuming pure hydrogen composition. For details on the HDUST code and the hydrodynamics of a VDD; see Carciofi & Bjorkman (2006, 2008) and Carciofi (2011).

3.1. Central star

The codeHDUST takes into account both the geometrical deformation of the star due to fast rotation and the latitude-dependent flux from the stellar surface due to gravity darkening. The stellar surface itself is divided into a number of latitude bins each with its own local gravity and temperature with a corresponding Kurucz (Kurucz 1979) synthetic spectrum assigned to it.

A rotating star is fully described by five parameters: mass M, polar radius R_p, luminosity L, rotation rate W, and gravity darkening exponent β. The star is a geometrically deformed isopotential surface in the Roche approximation based on the rotation rate W, defined as W = v_rot/v_orb, where v_rot is the rotational velocity at the stellar equator and v_orb is the Keplerian velocity of orbiting material just above the stellar equator (Sect. 2.3.1 of Rivinius et al. 2013a). In the case of critical rotation, the parameter W becomes unity. The gravity darkening law is adopted in the classical form T_eff ∝ g^β (von Zeipel 1924).

Two additional geometrical parameters are needed to be able to compare the model with observations: The distance d and the inclination angle i of the spin axis of the star with respect to the observer.

3.2. The disk

The VDD model in which turbulent viscosity is responsible for the disk growth and transport of the angular momentum outward has so far been able to explain most observational aspects of the Be star disks (Bjorkman 2012).

The VDD is flaring with the following dependence of the scale height H on the distance r from the stellar surface, if one assumes hydrostatic equilibrium and an isothermal disk $\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{H}}_{\mathrm{0}}{\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{e}}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}$(1)where R_e is the equatorial radius of the rotationally deformed star and ${\mathit{H}}_{\mathrm{0}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}{\mathit{v}}_{\mathrm{orb}}^{-1}{\mathit{R}}_{\mathrm{e}}$ is the scale height at the base of the disk. c_s = [(k_BT_k)/(μm_H)] ^1/2 is the isothermal sound speed, k_B is the Boltzmann constant, T_k is the gas kinetic temperature, μ is the mean molecular weight, and m_H is the mass of hydrogen atom. Haubois et al. (2012) studied the hydrodynamics of Be disks and explored how the disk surface density evolves with time in response to varying disk feeding rates. Even though a dynamically active disk can have a quite complex structure, a much simpler situation arises when the disk is fed at a constant rate for a sufficiently long time. In this case, the surface density structure of a VDD under the assumption of an isothermal gas and purely Keplerian circular orbits can be described by $\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{Ṁ}{\mathit{v}}_{\mathrm{orb}}{\mathit{R}}_{\mathrm{e}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{3}\mathit{\pi \alpha }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{\left[}\mathrm{\right(}{\mathit{R}}_{\mathrm{0}}\mathit{/}\mathit{r}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{]}}\mathit{,}$(2)where α is the viscosity parameter (assumed to be constant throughout the disk), Ṁ is the mass decretion rate, r is the distance from the stellar equator, and R₀ is an integration constant connected with the physical size of the disk (Bjorkman & Carciofi 2005). For r ≪ R₀, the dependence of Σ(r) becomes a power-law Σ(r) ∝ r^-2. The disk is in hydrostatic equilibrium in the vertical direction, which in the isothermal case corresponds to Gaussian distribution. The falloff of the volume density is then$\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\sqrt{\mathrm{2}\mathit{\pi }}\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\exp \left[\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{z}}{\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\right)}^{\mathrm{2}}\right]\mathit{,}$(3)therefore $\mathit{\rho }\mathrm{\propto }\mathrm{\Sigma }\mathit{/}\mathit{H}\mathrm{\propto }{\mathit{r}}^{\mathrm{-}\mathit{n}}\mathit{,}$(4)with the radial density exponent n equal to 3.5, however, the disk is not radially isothermal. Initially the temperature falls with radius to a minimum that occurs in the dense inner parts of the disk. The location of the minimum depends on the disk density; for higher densities it shifts outward. The temperature minimum results in a local change of the radial density falloff exponent n and of the flaring coefficient (Carciofi & Bjorkman 2008).

Numerical solutions (e.g., Okazaki et al. 2002) show that the disk is separated into two parts: a subsonic inner part (radial velocity v_r ≪ c_s), where the outflow is driven by viscosity; and a transonic outer part (v_r ≳ c_s), where the outflow is driven by the gas pressure. The photoevaporation radius where the transition occurs is the so-called transonic critical radius R_c, which is on the order of hundreds of stellar radii (Krtička et al. 2011). However, what is considered to be the outer radius of a Be disk depends on whether the star is in a binary system or not. If the star is in a binary system, the disk may be truncated at the tidal radius, which can be a complicated function of the binary parameters and of the viscosity parameter α (Panoglou et al. 2015). For isolated stars or well-detached binaries, however, we expect R_out ≈ R_c. Radio observations have to be used to investigate these outer parts of the disk.

3.3. Modeling procedure

For the present study, we have adopted two VDD models for the disk. First we start with a simplified analytical power-law model followed by the fully self-consistent solution of steady viscous decretion. As the latter model modifies the prescribed density structure according to the temperature solution, comparing the results of the two models allows us to test how the disk nonisothermality affects the observables. Improvement of the global fit to the observations after introducing the full VDD solution would then represent an indirect observational proof of the structural changes due to the disk nonisothermality.

The first model, which we dub parametric model, is based on the isothermal VDD density structure for which the vertical distribution is Gaussian and the radial distribution is a power law in the form ρ ∝ ρ₀r^{− n}, where ρ₀ is the density at the base of the disk. Although this corresponds to the theoretical density structure of an isothermal disk (when n = 3.5), the disk itself is not isothermal as its temperature is calculated in radiative equilibrium by HDUST. Values of ρ₀ were explored in the interval 1 × 10^-12–1 ×10^-10 g cm^-3, which is typical for Be disks. The scale height at the base of the disk H₀ is an additional parameter, which is set by the sound speed c_s and the disk orbital speed v_orb (see Eq. (1)). As v_orb is given by M, R_p, and W, only a proper choice of the kinetic temperature T_k (to calculate c_s) is needed to determine the scale height H₀. We chose the mass-averaged temperature of the disk to represent T_k.

Varying the density falloff exponent n allows us to explore a standard, isolated disk with isothermal density structure (n = 3.5) as well as a disk truncated by a binary companion, in which the parts inward of the truncation radius may attain a shallower density profile (n< 3.5). Density profiles associated with a disk that has not yet reached a quasi-steady-state configuration (n> 3.5, see Porter 1999; Haubois et al. 2012) were also explored (although this case is unlikely for β CMi, as its disk has been present and stable for a long period already). In our modeling, n was allowed to vary in the interval 3.0–4.0.

The second, self-consistent model calculates the full solution of steady viscous decretion following Carciofi & Bjorkman (2008). In this case, the radial and vertical density structures are solved to become consistent with the temperature solution and therefore the parameters n and H₀ are no longer needed. Instead of ρ₀ describing the base density, the ratio of the mass decretion rate and the viscosity parameter Ṁ/α needs to be specified (see Eq. (2)). The Ṁ/α parameter was chosen to correspond to the ρ₀ of the parametric model so that the resulting mass of the disk would remain approximately constant, thus allowing for a direct comparison between the models. The observables we study are only sensitive to the density scale of the disk, and do not carry any information about its radial velocity because it is much smaller than the sound speed. For this reason, the mass-loss rate Ṁ cannot be determined uniquely; we can only determine ρ₀ or equivalently the ratio Ṁ/α. Finally, there is an additional parameter in the form of the integration constant R₀, which was assumed to correspond roughly to the critical radius R_c given by the approximate relation of Krtička et al. (2011), i.e., $\frac{{\mathit{R}}_{\mathrm{c}}}{{\mathit{R}}_{\mathrm{e}}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{10}}{\left(\frac{{\mathit{v}}_{\mathrm{orb}}}{{\mathit{c}}_{\mathrm{s}}}\right)}^{\mathrm{2}}\mathrm{·}$(5)The parameter R_c, calculated from the best-fit model values, is ~450 R_e.

In both models, we also need to specify the disk outer radius R_out, which in the modeling context represents the outer edge of the grid used for the Monte Carlo simulation, beyond which matter is absent. This parameter allows us to explore possible truncation of the disk by a binary companion.

The basic characteristics and free parameters of the models are listed in Table 4.

3.4. Initial estimates for the central star parameters

β CMi has already been a subject of several detailed studies focusing on different observational and physical aspects. β CMi was previously studied by Saio et al. (2007), whose analysis of MOST satellite observations led to the first detection of nonradial pulsations in a late-type Be star with evidence for an almost critical rotation. Elsewhere, high-resolution spectrointerferometry and spectroscopy of β CMi was used to provide evidence of Keplerian rotation of the disk (Kraus et al. 2012; Wheelwright et al. 2012). The central star parameters adopted in these studies were used as initial constraints for our modeling.

The value of M = 3.5 M_⊙ was adopted in accordance with individual literature estimates (Kraus et al. 2012; Saio et al. 2007) and stellar evolution models for the spectral type B8. Although M influences the scale height of the disk as well as the critical rotation rate, changing it over the interval 3–4 M_⊙ has little overall influence on the observables, and it was therefore kept fixed at the value of 3.5 M_⊙.

We explored the possible values of R_p and L in the intervals 2.4–3.6 R_⊙ and 160–240 L_⊙, respectively. The rotation rate W was assumed to be 0.7–1.0, corresponding to the observed rotation rates of Be stars (Sect. 3.1 of Rivinius et al. 2013a).

Recent interferometric observations as well as theoretical developments indicate that the gravity darkening law in the classical form T_eff ∝ g^β with β = 0.25 is valid only in the limit of slow rotators, while in rapidly rotating stars β can be significantly lower. The recent model of Espinosa Lara & Rieutord (2011) in which the barotropic assumption (pressure dependent only on density) is relaxed is in general agreement with six interferometric measurements of the β parameter in rapid rotators, as is shown in Fig. 13 of Domiciano de Souza et al. (2014). In this work we treat the gravity darkening exponent β as a function of W, according to the model of Espinosa Lara & Rieutord (2011). However, in their model, gravity darkening is in fact not described by a power law (see their Eq. (31)).

3.5. Geometrical parameters

The distance d to β CMi is 49.6 ± 0.5 pc as measured by the Hipparcos satellite (van Leeuwen 2007). As varying d in such a low uncertainty interval has little effect on the model observables, d was kept fixed at the value of 49.6 pc. The inclination angle i is clearly of an intermediate value with the most precise estimate from modeling of AMBER observations being 38.5 ± 1° (Kraus et al. 2012). We explored the possible values of i in the interval of 35°–50°.

4. Results

In this section the ability of the VDD model, in the two particular formalisms adopted (parametric vs. self-consistent), to reproduce individual observations is detailed.

The fitting procedure was executed in the following way: First, the almost purely photospheric UV spectrum was used to constrain the central star parameters R_p and L. Next, the SED from IR to radio wavelengths, where the disk significantly contributes to the observed flux, was used to determine the disk parameters ρ₀ and n (or Ṁ/α). Finally, the spectral slope of the linear polarimetry allowed us to determine the rotation rate W (and β), while we used the interferometric shape extracted from the AMBER observations to constrain the geometrical parameter i. The resulting best-fit parameters for both models are listed in Table 5 and the derived stellar parameters are listed in Table 6.

The remaining observations were used as consistency checks for the parameters determined in the way described above. The polarization level allows us to check the consistency of i, ρ₀, and n in the inner parts of the disk. The equivalent widths of the different hydrogen lines provide us with another way to check the density structure in the various parts of the disk where the emission from the respective lines originates. The final consistency checks are provided by the observed interferometric visibilities, which measure the sizes of the corresponding emitting regions, along with the interferometric phase shifts, which depend on both the velocities and densities of the corresponding parts of the disk.

Comparing the outcomes of the best-fit (parametric and self-consistent) VDD models then allows us to explore the nonisothermal effects on the observables caused by structural changes in the disk introduced with the nonisothermal density solution. Finally, we look at the evidence in our data for the disk truncation and the presence of an unseen binary companion by exploring the outer parts of the disk using the radio data and several unexpected nonhydrogen spectral features.

4.1. SED

4.1.1. UV spectrum

The observed UV spectrum was found to represent the photospheric spectrum of the central star (Fig. 2) almost purely. Therefore the fluxes in the UV region are only influenced by the central star and can be used to constrain its parameters. As M is kept fixed and β is a known function of W, the parameters that remain to be determined are R_p, L, and W. While varying L in the interval 160–240 L_⊙ shifts the flux in the whole UV interval mostly in magnitude, varying R_p has a much larger effect on T_pole and hence on the spectral slope. With W (and the corresponding β) fixed at the value determined from the analysis of visual polarimetry (Sect. 4.3), and i fixed at the value determined from fitting AMBER interferometry (Sect. 4.4), R_p and L can be determined by fitting the UV spectral slope and the UV flux level, respectively. Out of the grid of models computed, the resulting best-fit parameters are R_p = 2.8 R_⊙ and L = 185 L_⊙. The model reproduction of the SED structure in the 1000–10 000 Å interval is shown in Fig. 3 (as the disk contribution is less than 3% in these parts on average, the outcomes of the two VDD model formalisms are identical).

4.1.2. Optical and IR SED

As discussed in the Introduction, the size of the disk pseudophotosphere grows with wavelength. Consequently, the near-IR regions are useful to probe the inner parts of the disk and constrain the disk base density ρ₀, while the SED structure at longer wavelengths traces the density falloff exponent n. The resulting ρ₀ from the best-fit parametric model is 2.0 × 10^-12 g·cm^-3; the implications for n in the light of the whole SED structure are discussed below.

The parametric model with a single value of n throughout the whole disk does not seem to provide an entirely satisfactory description of the density profile. In Fig. 4 we see that the parametric model with n = 3.5 reproduces the SED well only until ~15 μm, while it underestimates the far-IR and radio fluxes by about 40–70%. A shallower density profile with n = 3.0 (and the same ρ₀) reproduces the SED longward of ~15 μm much better at the expense of the mid-IR fluxes, which become slightly overestimated by about 20%, as may also be seen in Fig 5, which shows the comparison of the parametric model to the N-band MIDI photometry.

The self-consistent model SED is very similar to the parametric n = 3.0 model at wavelengths up to ~15 μm, showing only slight effects of nonisothermality in the inner, central parts of the disk, which are consistent with the predictions: The temperature initially decreases, which leads to a shallower radial density profile in the inner parts of the disk (see Fig. 4 of Carciofi & Bjorkman 2008). At longer wavelengths, the self-consistent model approaches the n = 3.5 model, as the disk temperature rises and becomes roughly isothermal in its more extended regions.

Overall, we find that neither the parametric model with a single n, nor the self-consistent model, are able to perfectly reproduce the observed SED structure, which suggests that the disk has a more complicated density profile than a simple power law, or steady decretion. We discuss the density structure further when we compare the model to the remaining observables, and we return to the issue of nonisothermality in Sect. 4.7.

4.2. The hydrogen emission line profiles

The overall shape, equivalent width, FWHM, and the ratio of maximum to continuum intensity of the observed Hα line profile is in reasonable (with respect to the used methods) agreement with the historic measurements from 1950 to 1991 (Slettebak & Reynolds 1978; Slettebak 1982; Andrillat & Fehrenbach 1982; Hanuschik et al. 1996). The hydrogen lines Hα, Hβ, Hγ, and Brγ all show evidence of a double-peaked emission component. Using the best-fit parametric and self-consistent models, we calculated both line profiles and images to provide an additional consistency check of the disk density structure, as the higher Balmer lines originate progressively closer to the star. The line profile results (Fig. 6) clearly show that a model with a single power law describing the radial density falloff is unable to simultaneously fit all four lines. While Hβ and Hγ, which originate in the inner regions of the disk, are reproduced well with n = 3.5, Hα and Brγ, which come from much more extended parts of the disk, require n = 3.0 to attain the observed level.

The outcome is similar to what was suggested by the SED analysis. A steeper falloff (higher n) is needed in the inner parts, where Hβ, Hγ, and near-IR continuum (but also the visual polarization, see the following section) originate, while in the more extended parts, where Hα, Brγ, and far-IR emission come from, a shallower density profile is needed to increase the size of the optically thick emitting area and hence the line equivalent width. This again indicates that a single power law throughout the disk is not a suitable description of the density profile. The outcome of the self-consistent model for these four lines is almost identical to the parametric model with n = 3.5.

For very optically thick emission lines, such as Hα, the line photons are absorbed and re-emitted many times within the Sobolev zone before they escape. This diffusion process greatly increases the probability that some fraction of the Hα line photons are scattered by free electrons, which have much larger thermal velocities than the hydrogen lines. To simulate this effect, we thermally broadened a fraction, f = 0.6, of the Hα photons by convolving those photons with a Gaussian of width v_es = 300 km s^-1. The remaining fraction (1−f) of the Hα photons are left unbroadened. The values of f and v_es were chosen empirically to best reproduce the wings of the observed Hα profile. The electron thermal velocity at T_k = 6500 K is 310 km s^-1, which is roughly what we determined empirically for v_es. The remaining Balmer lines (Hβ, Hγ and Brγ) are much weaker, so we applied no electron scattering broadening to those lines.

4.3. The optical polarization

The HPOL measurements show a level of variability, which is hard to quantify due to their large errors and scarce temporal coverage throughout the years 1990–2005. They are consistent though with a very low level of polarization in the range of 0.01–0.10%, suggesting no large-scale mass injection or dissipation episodes have occured over the past few decades. The OPD data covering approximately the last four years show a polarization level of ~0.02–0.07% in the BVRI Johnson filters, with less scatter compared to the HPOL data and also very little temporal variation (Table A.1).

For the modeling purposes, we used the averaged values of the OPD measurements. Despite the uncertainties being quite large, the data clearly show a positive spectral slope of the polarization level in the Paschen continuum (see Fig. 7). This is an unusual feature and is not seen in any of the theoretical polarization spectra shown by Haubois et al. (2014); in that work, low-density models have either flat spectra, indicative of the prevalence of Thomson scattering in the total opacity (i.e., low-density models), or negative slopes for the cases where H opacity contributes significantly to the total opacity (i.e., high-density models). In view of that, we explored which model parameter could be responsible for the positive slope. We used only the Paschen continuum polarization levels to compare the model and observed slopes to avoid the influence of Balmer and Paschen jumps.

Increasing the rotation rate W closer to the critical value was found to be the major factor causing the polarimetric slope to attain the unusual positive slope shown by the observations. The dependence of model polarization on W is shown in Fig. 8. Only rotation rates higher than about 0.90 produce a positive polarimetric slope, while rotation rates higher than W = 0.98 are required to reproduce the observations within the error bars. The dependence of the polarimetric slope on the rotation can be probably explained by the fact that the flux constrast between the pole and the equator gets stronger at shorter wavelengths (e.g., Fig. 6 of Rivinius et al. 2013a). Since polarization is the ratio between the scattered to direct startlight, it follows that polarization should be smaller where the flux constrast is the strongest. This effect is only valid for tenuous disks whose dominant opacity (Thomson scattering) is gray. Denser disks attain a negative slope reflecting the strong contribution of H to the total opacity. The finding of such a fast rotation is supported by the results of Saio et al. (2007), who found nonradial pulsation modes (prograde sectoral g-modes of m = −1) excited in β CMi, which they only predict to be stable if the star rotates nearly critically.

The predictions of the adopted models with W = 0.99, compared to the average of OPD measurements, are shown in Fig. 7. The parametric model with n = 3.0 overestimates the polarization level, while using n = 3.5 matches the observed level well. This further corroborates what was suggested by the SED and line profile analysis, i.e., a steep density falloff in the inner parts and a shallower profile farther on.

The polarization level is strongly dependent on the inclination angle i with the maximum level observed when i ~ 70−80° (Halonen & Jones 2013). However, the inclination was found to be well constrained by the AMBER data (see Sect. 4.4), which allows us to use the polarization level to constrain the density structure of the inner disk without a degeneracy with respect to i.

4.4. AMBER spectrointerferometry

AMBER data allow us to study the geometry and extension of the Brγ emitting region. To accomplish this, we calculated the synthetic Brγ images (several shown in Fig. 9) across the emission lines, depicting the corresponding isovelocity regions. From these images we determined the predicted AMBER visibilities and phase shifts across the Brγ line. AMBER’s high spectral resolution reveals the spectrointerferometric signatures in the differential visibilities and phases that are expected for circumstellar disks. As was reported by Kraus et al. (2012), and discussed theoretically by Faes et al. (2013), phase reversals occur in the differential phase signatures of β CMi across the Brγ line. Using the terminology of Faes et al. (2013), the disk of β CMi was observed outside the astrometric regime, which means that the ratio between the baseline length and the distance to the star exceeds 1.5 m pc^-1. Thus the differential phases are no longer proportional to photocenter displacement and phase signatures are more complicated than simple S-shaped profiles usually observed for marginally resolved disks.

Similar to the observables believed to originate in the more extended parts of the disk, a shallower density profile with n = 3.0 is needed to reproduce the AMBER observations, while the predictions of the n = 3.5 parametric model and the self-consistent model are almost identical, clearly underestimating the size of the line emitting region (Fig. 10). The fit of the n = 3.0 parametric model to the two sets of AMBER data is remarkable, with a reduced χ² = 1.20. This kind of good agreement again implies that the disk velocities are Keplerian, which is strong support for the VDD model. The inclination is well constrained by the shape implied by the visibilitiy data along the various baseline orientations with the best-fit value being $\mathit{i}\mathrm{=}\mathrm{4}{\mathrm{3}}_{\mathrm{-}{\mathrm{2}}^{\mathrm{◦}}}^{\mathrm{+}{\mathrm{3}}^{\mathrm{◦}}}$. The position angle, PA, of the disk major axis projected on the sky was found to be 133$\begin{array}{c}\mathrm{+}\mathrm{4}{\mathrm{◦}}^{}\\ \mathrm{-}\mathrm{3}{\mathrm{◦}}^{}\end{array}$ (from north to east) in a reasonable agreement with the value of 140.0 ± 1.7° found by Kraus et al. (2012). The average polarization position angle from the HPOL synthetic UBVRI filter data and OPD measurements is 41.6±1.8° and 50.0±17.8°, respectively, in a good agreement with the interferometric PA, as the polarization direction is usually perpendicular to the disk major axis, but small variations can be expected from the disk asymmetries responsible for the V/R variations.

4.5. NPOI interferometry

The NPOI data, which contain the Hα line emission in its 15 nm-wide channel, provide an important consistency check of the modeling of the Hα line profile, which could only be reasonably reproduced with a shallow, n = 3.0, radial density falloff. Indeed, the NPOI visibilities are better reproduced with the n = 3.0 model (χ² = 1.40), revealing a larger Hα emitting region than expected from the n = 3.5 model (χ² = 2.08). The predictions of the two models are compared with the data in Fig. 11. We used an average angular radius of the central star to calibrate all the PAs, while the proper approach would be to use the exactly corresponding angular radii to calibrate the individual PAs. However, we expect that this would bring only minor improvement of the overall good fit.

4.6. FSU-A, MIDI, PIONIER, and CHARA interferometry

The remaining interferometric measurements are either very scarce in their uv-plane coverage or too noisy to provide a reliable consistency checks of the best-fit model parameters. Nevertheless, the model reproduces these data reasonably well (Table 7).

4.6.1. VLTI

The fit to the three FSU-A K-band visibilities is shown in Fig. A.1. The MIDI data covering the N-band spectral region are consistent with an unresolved disk, i.e., values of visibilities ~1 within 1–2σ in a good agreement with both the n = 3.5 and n = 3.0 parametric models (figure not shown). The only usable outcome of the PIONIER observations was a closure phase measurement consistent with zero (± 1.5°), which eliminates any companion up to 5% flux contrast in the 1.5–50 mas separation range.

4.6.2. CHARA

Despite the MIRC data having the best uv-plane coverage and the longest projected baselines of all the interferometric measurements presented, significant scatter of the data, often of measurements taken during the same night for almost identical baselines prevented the data from being of much use (Fig. 12). The cloud of residuals is not centered on zero, but below, i.e., the observations often show a lower fringe contrast than predicted. Diminished fringe contrast could be explained by observing issues and/or weather conditions that were not accounted for.

Finally, the model comparison to the K’-band CLIMB data is plotted in Fig. A.2.

4.7. Nonisothermal effects on the observables

So far, we have mostly focused on presenting the results of the parametric model as it allows us to test alternative density profiles. In this section, we compare the outcomes of the parametric model with n = 3.5 and the self-consistent model and discuss the nonisothermal effects on both the disk structure and model observables that arise when the density structure is self-consistently calculated together with the temperature structure.

The temperature, level populations, and density structures of the two models along the disk midplane are plotted in Fig. 13. The temperature is found to decrease up to about 3 R_e, then to steeply rise, and from about 5 R_e further the disk becomes isothermal and almost neutral. Introducing the self-consistent solution leads to a slight compression of the disk in the vertical direction, while the radial density structure changes as expected from the temperature structure, i.e., the falloff is shallower where the temperature gradient is negative and steeper where the temperature gradient is positive (Carciofi & Bjorkman 2008). This has an effect on the near-IR SED structure, as shown in Fig. 4, and also on the polarization level, which slightly increases because of enhanced density in the inner parts.

The effect of the self-consistent solution on the density structure and SED is opposite to what was inferred from the SED and spectral line analysis, a steep density falloff followed by a shallower falloff. This appears to indicate that some mechanism, other than nonisothermality, is modifying the disk structure. Note, however, that the non-isothermal effects on the observables are small, as expected for a low-density disk around a late-type B star (Carciofi & Bjorkman 2008). The presence of a possible unseen binary companion would most likely have stronger effects on the disk structure, and would prevent us from finding a global improvement of the model fit to the observations when making the model fully self-consistent.

It is also interesting to note that the model Brγ visibilities and phases are identical for the parametric model with n = 3.5 and the self-consistent model, which means that while the density structure changes locally, the outcome for the whole Brγ-emitting region for these two models is the same. We next turn to the evidence in our data for the presence of an unseen binary companion.

4.8. Estimating the disk size from radio fluxes

Contrary to the near-IR, which probes the inner disk, the radio fluxes probe the outermost reaches of the disk. For an isolated disk, the disk can grow to the photoevaporation radius. However, if the disk is physically truncated by a binary, the pseudophotosphere effective radius eventually reaches the disk physical size and the disk flux excess no longer increases with wavelength. Consequently, the SED presents an inflection at the wavelength where ${\overline{\mathit{R}}}_{\mathit{\lambda }}\mathrm{\approx }{\mathit{R}}_{\mathrm{out}}$, and follows the photospheric spectral slope at longer wavelengths (see Fig. 11 of Vieira et al. 2015, for an illustration). Therefore, the observation of this kind of transition in the SED slope can place important constraints on the disk physical size.

Under the assumption that the shallower density profile needed to reproduce the mid- to far-IR SED structure also holds throughout the radio emitting region, we tested which disk radii reproduce best the mm and cm SED structure. In the upper panel of Fig. 14, we plot the sub-mm to cm SED for models with different R_out and, in the lower panel, the corresponding reduced χ². The APEX and JCMT fluxes are already reproduced for R_out = 20 R_e and increasing R_out beyond this value does not lead to further enhancement of the model flux. The physical explanation for this is that the size of the pseudophotosphere at these wavelengths is about 10–20 R_e, so increasing the disk size beyond this value does not result in an increase in the disk emission. The behavior at CARMA wavelength is similar, but R_out of at least 40 R_e is needed. At 2 cm, however, the situation is different. While the lowest residual between the model and observed flux is seen for R_out = 35 R_e, a higher value of R_out overestimates the observed flux. The best-fit value of R_out is found to be 35$\begin{array}{c}\mathrm{+}\mathrm{10}\\ -5\end{array}$R_e.

In Fig. 15 we show the synthetic images and graphical representation of the brightness profiles at wavelengths corresponding to the APEX, CARMA, and VLA observations with the estimates of the corresponding pseudophotosphere sizes. Indeed, as expected from the model behavior, at the shorter wavelengths ${\overline{\mathit{R}}}_{\mathit{\lambda }}$ is between 10–20 R_e. However, ${\overline{\mathit{R}}}_{\mathit{\lambda }}$ at 2 cm is about 40 R_e, thus larger than the estimated disk size. The 2 cm results can therefore be interpreted as a clear sign of truncation by an unseen binary companion, however, this conclusion rests on the reliability of a single non-contemporaneous data point.

4.9. Spectral features of other elements

To explore the additional spectral features in the interval covered by the FEROS and ESPaDOnS spectra, we used an updated version of the BRUCE3 code (see Rivinius et al. 2013b, and the references therein) to compute a purely photospheric, high-resolution spectrum for the stellar parameters listed in Table 5.

An unusual feature is present in both the FEROS and ESPaDOnS spectra. The calcium triplet (CaII λ8498, 8542, and 8662 Å) is in emission (blended with the Paschen lines Pa16, Pa15, and Pa13), but is noticeably stronger than the emission from the unblended Paschen line Pa14 (Fig. 16). This effect was linked to possible binarity of Be stars by Polidan (1976), while a correlation of the presence of this feature in a number of B subdwarf (sdB) stars with a cool companion was reported by Jeffery & Pollacco (1998). The emission in this case probably originates in the very outer disk, which is possibly partially accreting onto the secondary component.

Moreover, in the CaII λ3934 line, the emission component is noticeably narrower than the absorption component (left panel of Fig. 17), while the opposite is expected for a line originating in the fast-rotating inner disk. Therefore, similarly to the calcium triplet mentioned above, the emission seems to actually originate in the outer disk, close to a possible secondary companion.

Finally, the ultraviolet CIV λ1548 line exhibits a P Cygni profile with a blue edge velocity of ~230 km s^-1 (right panel of Fig. 17), which is a clear signature of a fast stellar wind. The presence of this line is not expected in this late-type B star (see Fig. 11 of Rivinius et al. 2013a). Additionally, the observed velocity is lower than is typical for late B stars. This suggests that this kind of wind might be excited by the radiation of a companion, which interacts with the outer disk of the primary Be star.

5. Discussion on the possible binarity

Using the radio flux measurement at the wavelength of 2 cm in combination with a well-sampled SED at shorter wavelengths, we found observational evidence suggesting truncation of the disk of β CMi at about 35 R_e. The most probable mechanism that can truncate the disk at this radius is the tidal influence of an orbiting binary companion. The presence of an unseen binary companion is also suggested by weak V/R oscillations of the Hα line and the appearance of the CaII λ3934 emission line and the triplet CaII λ8498, 8542, 8662. Additionally, the CIV λ1548 line sugggests a hot, possibly subdwarf (sdB or sdO) companion. The proposed companion has to be faint, as its signatures are seen neither in the spectra nor in the SED. The detection limits from interferometric measurements set the spectral type of a main-sequence companion to be F5 or later, thus putting an upper limit of ~0.5 on the mass ratio q.

To test whether a binary companion could truncate the disk at the radius estimated from the modeling, we made the assumption that the V/R variations in Hα are indeed caused by a companion, in which case the period of these oscillations (~183 days, Folsom et al. 2015) can be used as an estimate of the orbital period. We then performed simulations of the binary system in a circular orbit with a smoothed particles hydrodynamics (SPH) code (Okazaki et al. 2002), which assumes an isothermal disk, to determine the truncation radius, R_t, resulting from the presence of the possible companion. In the SPH simulation, the Be star disk is initially nonexistent. We then assume a constant mass injection rate from the equator of the star onto the equatorial plane, so that after letting the system evolve, we achieve a base density ρ₀ equal to its constrained value (Table 5). We explored the outcome of those simulations (shown in Fig. 18) for different viscosity parameters, α ∈ { 0.1,1.0 }, and mass ratios, q ∈ { 0.1,0.5 }. Figure 18 shows the surface density structure of representative models after a long period of disk evolution (30 orbital periods) to ensure that the disk has reached a quasi-steady-state configuration (Panoglou et al. 2015; Haubois et al. 2012). Following Okazaki et al. (2002), the truncation radius is estimated by fitting the expression $\mathrm{\Sigma }\mathrm{=}\mathit{A}\frac{\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{t}}{\mathrm{)}}^{\mathrm{-}\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{t}}{\mathrm{)}}^{\mathit{l}}}$(6)to the results of the SPH simulations, where A, R_t, k, l are constants to be fit. The resulting truncation radii range from ~18 R_e (α = 0.1; q = 0.5) to ~25 R_e (α = 0.1; q = 0.1), as compared to an R_out of 35$\begin{array}{c}\mathrm{+}\mathrm{10}\\ -5\end{array}$R_e determined from the SED fitting. A certain discrepancy is to be expected because in the HDUST model, the disk is entirely cut off at ~35 R_e. In contrast, the truncation radius in the SPH simulations represents a point were the density falloff becomes steeper (Fig. 18), but the disk still exists beyond this point. This means that the size of the pseudophotosphere at 2 cm should in principle be somewhat larger than R_t, which would reconcile the two results.

Out of the explored mass ratios and viscosity parameters, the q = 0.5 and α = 0.1 simulation is the only one that clearly causes the density structure of the disk inward of R_t to attain a shallower profile of n = 3.0, which is the value that was determined for intermediate part of the disk (inward of R_out) from the HDUST modeling. However, the steeper density profile in the very inner parts of the disk is not reproduced by the simulation, although there is a possibility that the introduction of nonisothermality into the SPH simulation might show the desired effect.

In any case, although the proposed binarity is quite probable, this must be viewed with caution. Firstly, the evidence for the disk truncation is based only on one single measurement from more than 30 years ago. Despite the apparent stability of the disk over the last ~15 years, as shown by our spectroscopic data set, and an overall stability over the last ~65 years according to the literature, we cannot exclude the possibility that the mass injection rate into the disk has changed over the last decades to a degree that would destroy evidence of truncation. Secondly, it is unclear whether the complicated density structure including the inner parts is compatible with the presence of the binary companion, although the tidal interaction seems to offer the most straightforward explanation. A careful analysis of nonisothermal SPH simulations is necessary to explore if the presence of the companion offers the explanation, or if other causes are needed. Since this is beyond the scope of the present analysis, we plan to return to this issue in a follow-up paper. In that paper, we intend to study the density profile in more detail as well as execute our own analysis of the V/R oscillations using all available measurements of Hα and possibly other emission lines.

In addition to binarity, possible mechanisms for truncating the disk are photoevaporation by gas pressure and line driven radiative ablation (Feldmeier et al. 1999; Gayley et al. 1999). The former was already examined in Sect. 3.2 and affects the disk structure only at very large radii. The latter, however, is a possibility to consider. Using the stellar wind mass-loss rate estimate of Castor et al. (1975), Krtička et al. (2011) have determined that the disk wind mass-loss rate scales with the stellar wind mass-loss rate, but with values lower by more than an order of magnitude. This implies that the disk ablation produces a mass loss much lower than that of a main-sequence star. As the winds of late-type main-sequence B stars are typically weak, we conclude that for stars, such as β CMi, wind ablation should not play a role in the truncation of the disk.

6. Conclusions

We tested the ability of the VDD model to reproduce arguably the largest data set of multiwavelength and multitechnique observations of a Be star (β CMi), including the SED covering five orders of magnitude in wavelength, emission line profiles of four hydrogen lines, visual polarimetry, high-resolution, near-IR spectrointerferometry and visual, near-IR, and mid-IR broadband interferometry. Neither a single power-law structure nor steady viscous decretion appears to be an entirely adequate description of the density profile. Observables originating in the inner disk (polarization, higher Balmer lines, near-IR photometry) require the steeper radial power-law exponent n = 3.5, while observables originating from more extended regions (Hα, Brγ, mid- to far- IR photometry, sub-mm/mm photometry) require the shallower n = 3.0. It is shown that nonisothermal effects on the disk structure cannot account for this more complicated density behavior. Instead, an outer tidal influence most probably exerted by an unseen binary companion is needed, which remains a possibility to be investigated in detail in the future.

Analysis of radio observations suggests truncation of the disk at the distance of 35$\begin{array}{c}\mathrm{+}\mathrm{10}\\ -5\end{array}$R_e from the star. This value was compared with the truncation radius exerted by a possible binary companion with the orbital period corresponding to the period of weak V/R oscillations of the Hα line with an overall good agreement. Together with spectroscopic evidence and a nondetection by any of the other observations, either a late (F5 or later) main-sequence companion or a subdwarf (sdB or sdO) companion is favored.

We report on the diagnostic potential of the polarimetric slope in the Paschen continuum for constraining the rotation rate W. The observed positive slope could only be reproduced with rotation rates very close to critical, with a favored value of W ≳ 0.98. The almost purely photospheric UV spectrum allowed us to constrain the polar radius R_p = 2.8 R_⊙ and luminosity L = 185 L_⊙ of the central star. Finally, AMBER spectrointerferometry facilitated a precise estimation of the inclination angle ($\mathit{i}\mathrm{=}\mathrm{4}{\mathrm{3}}_{\mathrm{-}{\mathrm{2}}^{\mathrm{◦}}}^{\mathrm{+}{\mathrm{3}}^{\mathrm{◦}}}$). Therefore, to conclude, these results underscore the utility of multitechnique and multiwavelength observations to constrain the disk structure unambiguously, as well as the use of radio observations to study the properties of the outer parts of Be disks.

Online material

Appendix A: Additional material

---

Acknowledgments

R.K. would like to acknowledge the kind help of Dietrich Baade in the final stages of preparation of the paper. R.K. thanks Petr Harmanec for being the observer of several spectra used in this work. The research of R.K. was supported by grant project number 1808214 of the Charles University Grant Agency (GA UK) and by the research program MSM0021620860 (MŠMT ČR). A.C.C. acknowledges support from CNPq (grant 308985/2009-5). D.M.F. acknowledges support from FAPESP (grant 2010/19060-5). R.G.V. acknowledges the support from FAPESP (grant 2012/20364-4). D.P. acknowledges support from FAPESP grant No 2013/16801-2. M.C. acknowledges, with thanks, the support of FONDECYT project 1130173 and Centro de Astrofísica de Valparaíso. This work has made use of the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was made possible by the Brazilian agency FAPESP (grant 2009/54006-4) and the INCT-A. Support for CARMA construction was derived from the states of California, Illinois, and Maryland, the James S. McDonnell Foundation, the Gordon and Betty Moore Foundation, the Kenneth T. and Eileen L. Norris Foundation, the University of Chicago, the Associates of the California Institute of Technology, and the National Science Foundation. Ongoing CARMA development and operations are supported by the National Science Foundation under a cooperative agreement, and by the CARMA partner universities. The Navy Precision Optical Interferometer is a joint project of the Naval Research Laboratory and the US Naval Observatory, in cooperation with Lowell Observatory and is funded by the Office of Naval Research and the Oceanographer of the Navy. C.T. would also like to thank student assistants, Bryan Demapan, Brennan Kerkstra, and Sandeep Chiluka for assistance with NPOI data reductions. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts.

References

All Tables

All Figures