© ESO, 2017

1. Introduction

Be stars offer unique possibilities for studying circumstellar disks. Because they are common among the bright, nearby stars – 17% of B-type stars are Be stars (Zorec & Briot 1997) – their disks are among typical targets for modern optical/near-infrared interferometers that have resolved them at the milliarcsec level. As a result, the structure of the inner parts of the disks is now well understood in the framework of the viscous decretion disk (VDD) model, first proposed by Lee et al. (1991) and further developed by, for example, Bjorkman (1997), Porter (1999), Okazaki (2001), and Bjorkman & Carciofi (2005). The central stars rotate close to break-up velocities, and an uncertain mechanism – the so-called Be phenomenon – acts in addition to rotation, leading to episodic or continuous mass ejection from the stellar equator. Subsequently, outflowing, ionised, purely gaseous disks, rotating in a nearly Keplerian way are formed (see Rivinius et al. 2013, for a recent review). In the VDD model, it is the turbulent viscosity that is responsible for the transport of the angular momentum outwards, and therefore for the growth of the disk.

In the last decades, observational techniques such as spectroscopy, linear polarimetry, and optical/near-IR interferometry have been used to constrain the physical structure of VDDs. Combining polarimetric and interferometric measurements, the disk-like structure of the circumstellar matter was unambiguosly confirmed (Quirrenbach et al. 1997) and the rotation law was subsequently confirmed to be nearly Keplerian (Meilland et al. 2007; Wheelwright et al. 2012; Kraus et al. 2012). However, it is important to note that these results apply to the inner parts of VDDs only (≲ 20 stellar radii, Vieira et al. 2015), that are responsible for the bulk of the disk emission in the optical and infrared (IR), as well as for the Balmer lines (Carciofi 2011). At larger radii, the disks are detectable only at radio wavelengths.

In this work, we are interested in the overall density structure of VDDs, focusing on the parts outwards of ~20 stellar radii. We constrain the density structure by compiling the observed spectral energy distribution (SED) from ultraviolet (UV) up to radio wavelengths. Due to the nature of the disk continuum emission – bound-free and free-free emission from the ionized hydrogen in the disk – fluxes at longer wavelengths originate from progressively larger areas of the disk. At the UV wavelengths, the SED usually consists of the stellar photospheric flux only, although in shell stars (Be stars seen close to edge-on), the UV photospheric spectrum can be dimmed by the surrounding disk and also strongly blanketed by the disk line absorption (mostly Fe). At longer wavelengths, the disk contribution to the SED becomes more significant until it dominates the SED in the mid-IR to radio domains. This is well explained in terms of the gas in the disk forming an optically thick pseudo-photosphere, whose size grows with increasing wavelength (see Vieira et al. 2015, for more details on the pseudo-photosphere concept applied to the continuum emission of Be disks).

Radio observations longwards of sub-mm wavelengths are of special interest for studying the outer parts of VDDs, as the disk pseudo-photospheres are larger at these wavelengths. In this work, we revisit the Be stars that were detected by the Very Large Array (VLA) observations at 2 cm by Taylor et al. (1990, six detected stars) and subsequently observed at mm wavelengths using the James Clerk Maxwell Telescope (JCMT) by Waters et al. (1991). The main outcome of the radio measurements was the discovery of an SED turndown, that is, a steepening of the spectral slope somewhere between mid-IR (10–60 μm, as measured by the IRAS satellite) and radio wavelengths, in the studied stars. Using a simple wedge-shaped disk model seen pole-on, Waters et al. (1991) were able to reproduce the observed SEDs with the disks having power-law density structure in the form ρ ∝ r^{− n}, with the exponent n ranging from 2.0 to 3.5. However, to reproduce the SED turndowns, disks truncated at a certain radius (ranging from 26 to 108 stellar radii) had to be used. Although the VDD model was not yet established at that time, it was already suspected that Be envelopes have disk-like shapes that originate from the central star and flow outwards. If the disks are dominated by rotation rather than outflow (as in the VDD model), it would follow that the most likely reason for the disk truncation is the presence of an (unseen) binary companion. However, because Waters et al. (1991) assumed wind-like outflows rather than disks dominated by rotation, they found truncation to be an unlikely explanation for the SED turndown, favoring other possibilities, such as changes in the disk geometry or additional acceleration of the disk material in the outer disk.

In the last decade or so, the VDD model has been firmly established as the framework in which most of the observable properties of Be disks can be explained (at least in their inner parts). Successful VDD reproductions of multi-technique and multi-wavelength observations of individual objects include, for example, ζ Tauri (Carciofi et al. 2009), β CMi (Klement et al. 2015, 2017) and 48 Lib (Silaj et al. 2016). This allows us to interpret the radio observations on a firm physical basis and therefore distinguish between the suggested scenarios to explain the SED turndown.

The goal of the present work is to carefully model the compiled SEDs of the Be stars detected in the radio by Taylor et al. (1990) using the VDD model implemented in the Monte Carlo radiative transfer code HDUST (Carciofi & Bjorkman 2006, 2008). For four of these stars, we have obtained new multi-band measurements in the cm domain, allowing us to study the properties of the SED turndown in detail. With the results of this study we address the question of the origin of the SED turndown, that is, whether it can be explained in the framework of the VDD model, which may be truncated by a binary companion, or if further physical ingredients in the model are needed.

In the following section we focus on the theoretical predictions concerning the structure of the outer disks of Be stars. In Sect. 3 we describe the observations used, and in Sect. 4 we explain the VDD model and the modeling procedure. The following section details the results for the six individual objects and the conclusions follow.

2. The outer parts of Be star disks

In an isolated VDD, the outflow velocity in the inner disk is highly subsonic and grows linearly with the distance from the central star (Okazaki 2001). At a certain point, sometimes called the critical or photo-evaporation radius, the outflow velocity reaches the local sound speed, which marks the transition between a subsonic inner part and a supersonic outer part of the disk. Outwards from this point, it is no longer viscosity but rather the gas pressure that drives the outflow, and the disk becomes angular momentum conserving. This transition may offer an explanation for the SED turndown, as the density decreases much more rapidly past the transonic part. However, the transition is only expected to occur at very large distances from the central star. An approximate relation for the critical radius R_c was derived by Krtička et al. (20111) for isothermal disks: $\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}}\mathit{,}$(1)where R_e is the stellar equatorial radius, v_orb is the disk orbital speed above the stellar equator and c_s is the speed of sound. The typical values of R_c are approximately 350 and 430 R_∗ for spectral types B0 and B9, respectively. According to the pseudo-photosphere model of Vieira et al. (2015), for typical disk densities, the pseudo-photosphere radius attains such large values only at wavelengths longer than 10 cm. Unless the disk is extremely dense, this means that only wavelengths longwards of 10 cm can probe the transonic regions. So far, no detections of Be disk emission at these wavelengths have been reported. It follows that the SED turndown is likely unrelated to the existence of a transonic regime in Be disks.

The only known physical mechanism that can truncate the disk at a closer distance to the central star is the tidal influence from an orbiting secondary companion. The effects of such a scenario on the VDD of the primary have recently been explored by smoothed particle hydrodynamics (SPH) simulations, under the assumption that the orbit of the secondary is aligned with the disk (Panoglou et al. 2016). For circular orbits, the influence of the secondary has two important effects on the Be disk. The first is the truncation of the disk at a distance close to the 3:1 resonance radius, which is much smaller than the orbital separation. Outwards of the truncation radius, the density starts to decrease much more rapidly, but these parts can still contribute to the emergent spectrum. The second is the so-called accumulation effect, which causes the parts of the disk inwards of the truncation radius to have a slower density fall-off (i.e., a smaller radial power law exponent). The accumulation effect has a larger impact with decreasing orbital period, decreasing viscosity, and increasing binary mass ratio. Considering the possibility that the orbit of the secondary is misaligned with respect to the equatorial plane of the disk, both the truncation and accumulation effects become much weaker (Cyr et al. 2017).

The incidence of binarity among Be stars remains an important open issue. The presence of a close binary companion was, in the past, believed to be the cause of the Be phenomenon, although that is not likely the case (Oudmaijer & Parr 2010). However, the evolution of a mass-transferring binary may produce a low-mass subdwarf O or B star (sdO/sdB, the mass donor that was originally the more massive component and is now stripped of its outer layers) and a fast-spinning star (the mass and angular momentum gainer, that is now the primary star). It is expected that at least some Be stars may have formed in this fashion. Owing to their lower luminosity, sdO/sdB companions are difficult to detect and only five Be+sdO/sdB systems have been reported so far: ϕ Per (Gies et al. 1998), FY CMa (Peters et al. 2008), 59 Cyg (Peters et al. 2013), o Pup (Koubský et al. 2012), and HR 2142 (Peters et al. 2016). They were revealed by periodic signatures in the UV spectra and by traveling emission components in the emission lines of the Be star caused by radiative interaction of the hot sdO/sdB star with the outer disk of the Be star. Other common companions of Be stars are neutron stars, which may emit X-rays as they accrete some of the matter from the disk of the Be primary. Other compact companions should include black holes and white dwarfs, although only one such system has been confirmed (Be star + black hole system HD 215227; Casares & Jonker 2014). The companions that are hardest to detect, though they are possibly common, are low-mass main-sequence stars.

Searching for binary companions using radial velocity (RV) shifts of spectral lines is complicated by their large rotational broadening. The exceptions from this rule are shell stars, which have narrow absorption lines in the centers of broad emission lines formed in the disk and thus orbiting companions are more easily detectable. Searches for binary companions of Be stars using speckle interferometry (Mason et al. 1997) and direct imaging with adaptive optics (Roberts et al. 2007; Oudmaijer & Parr 2010) led to the conclusion that the incidence of binaries among Be stars is around 30%. However, these results apply only for companions that are farther than 0.1 arcsec (~ 20 au for the most nearby targets), with a maximum magnitude difference of 10 mag in the I and K-bands. That means that Be companions such as main sequence objects of spectral classes F and cooler or sdO/sdB stars are not detectable with these techniques. Such binary companions would remain invisible, but could, in principle, be revealed through their influence on the density structure of the disk of the primary. If prevalent, the truncation of Be disks by such companions could be the cause of the observed turndown in the SEDs. Considering additional indirect evidence, line emission in the IR Ca triplet has been suggested to indicate binarity in Be stars (see, e.g., discussion by Koubský et al. 2011). Also, in known Be binaries in a close orbit (e.g., γ Cas), the Hα line often shows no clear peaks or symmetry, which is thought to indicate disturbances in the disk due to the orbiting companion.

3. Observations

The primary dataset used for the modeling is the SED measurements, namely the absolutely calibrated UV spectra and photometric measurements from visual to radio wavelengths. The secondary dataset contains visual spectra and polarimetry. The spectra were used to search for RV variations that could be assigned to the presence of an orbiting companion. The emission line profiles and the polarimetry were used to check for variability in the disk throughout the last decades.

The epochs, wavelengths, and references of the IR, sub-mm, and radio photometric measurements are given in Table 1. The epochs of the observations span over three decades. Since Be disks are often highly variable, this may introduce errors in our solutions. Nevertheless, the data for each target were still combined into a single dataset, which then represents average properties of the disk over the last decades. For more details on the variability of individual targets, and how it affects the solution, we refer to Sect. 5.

3.1. Ultraviolet spectra

For the UV part of the spectrum, data from the International Ultraviolet Explorer (IUE)¹ and the Wisconsin Ultraviolet Photo-Polarimeter Experiment (WUPPE)² were used. For the IUE spectra, large aperture measurements were used in all cases due to the problems with absolute flux calibration of the small aperture measurements. When available, we preferred low-dispersion spectra. The data show very little or no temporal variation for all stars and were therefore averaged for modeling purposes.

3.2. Visual and IR photometry

To compile all the available visual and IR photometric measurements, we made use of the Virtual Observatory SED Analyser³ (VOSA; Bayo et al. 2008) and the VO Spectral Analysis Tool⁴ (VOSpec; Arviset et al. 2008). These tools search the available catalogs and compile the SED for a given object. The visual domain catalogs that VOSA and VOSpec searched and that we subsequently used in this study include those of Ducati (2002), Hipparcos (Perryman et al. 1997), Mermilliod & Mermilliod (1994) and Tycho-2 (Høg et al. 2000). For the IR domain, data from the space mission MSX (MSX6C Infrared Point Source Catalog; Egan et al. 2003) was used.

In the near-IR domain we made use of JHKLMN magnitudes measured at three observatories: the University of Calgary Rothney Astrophysical Observatory (RAO) 1.5 m telescope, the National Optical Astronomy Observatory (NOAO), Kitt Peak 1.3 m telescope and the MKO 0.6 m telescope (Dougherty et al. 1991a). In the mid-IR domain, we used available data from IRAS (Helou & Walker 1988), AKARI (AKARI/IRC mid-IR all-sky Survey; Ishihara et al. 2010), and WISE (Cutri & et al. 2014) satellites. The catalog flux values were color-corrected in the same way as described by Vieira et al. (2017). For β CMi we also used data from the Spitzer Space Telescope (SST; Su et al. 2006). IR measurements corresponding to upper limits were discarded.

We note that for the star β Mon A, the visual and IR photometry actually corresponds to the combined photometry of the three components, all of which are Be stars of similar spectral types. This may introduce additional errors in our solution.

3.3. Radio observations

3.3.1. Sub-millimeter observations

For the star β Mon A, we used a new APEX measurement from the bolometer camera LABOCA (Siringo et al. 2009) at the wavelength of 870 μm. The observed flux is 23.5 ± 3.5 mJy. The separation between the A and B components is 7.1 arcsec. When centered on the A component, the B component falls inside the FWHM of the LABOCA beam, which is equal to 19.2 ± 3 arcsec. Although close to the edge of the FWHM, the B component may still slightly contribute to the observed flux.

For β CMi, we also have an APEX/LABOCA measurement, that was originally published in Klement et al. (2015). However, we newly reduced the data using Crush⁵ version 2.32-1 (Kovács 2008), with the updated flux value for β CMi being 38.6 ± 3.1 mJy. This value remains consistent within the error bars with the value used previously. For a more detailed description of the APEX/LABOCA measurements, see Sect. 2.2 of Klement et al. (2015). For β CMi, we also use the flux at 3 mm as measured by CARMA (see Table 1 of Klement et al. 2015).

3.3.2. Millimeter observations

At mm wavelengths, we used the previously published JCMT measurements of Waters et al. (1989, 1991) and the IRAM measurements of Wendker et al. (2000).

3.3.3. Centimeter observations

The previously published radio fluxes at cm wavelengths were taken from the VLA observations presented by Taylor et al. (1987, 1990), Apparao et al. (1990), Dougherty et al. (1991b) and Dougherty & Taylor (1992).

New radio observations of four stars were obtained with the Karl G. Jansky Very Large Array (JVLA) of the National Radio Astronomy Observatory in October, 2010 (Table 2). Each star was observed at four different wavelengths (6.0 cm, 3.5 cm, 1.3 cm, and 0.7 cm) in the C configuration (baselines ~ 0.035–3.4 km). For all four observing bands, the WIDAR correlator was configured using a setup that provided two sub-bands, each with a 128 MHz bandwidth, 64 spectral channels, and dual circular polarizations. Typical on-source integration times for each star were ~65, 13, 4, and 3.5 min at 6.0 cm, 3.5 cm, 1.3 cm, and 0.7 cm, respectively, with the exception of γ Cas, where a second observation of comparable duration was obtained at each of the three shortest wavelengths. Observations of each star were interspersed with observations of a neighboring bright point source to serve as a complex gain calibrator, and each session included an observation of 3C 48 at the appropriate observing band(s) to serve as an absolute flux density and bandpass calibrator. Data reduction was performed using the Astronomical Image Processing System (AIPS; Greisen 2003) and followed standard procedures for continuum data obtained with the WIDAR correlator as described in the AIPS Cookbook⁶. The data were loaded into AIPS directly from the archival science data model (ASDM) files using the Obit software package (Cotton 2008). However, the default calibration (“CL”) tables were recreated to update the gain and opacity information. The antenna positions were also updated to the best available values.

For all of the data obtained for this program it was necessary to perform an initial delay correction using a short (~1 min) segment of data from the flux/bandpass calibrator. Following the excision of poor-quality data, the data weights were calibrated based on the system power measurements from each antenna using the AIPS task tyapl. The absolute flux density scale was determined using the time-dependent 3C 48 flux density values from Perley & Butler (2013). Calibration of the bandpass and the frequency-independent portion of the complex gains was subsequently performed using standard techniques. At the two highest frequencies, the gain calibration was performed in two steps, first solving only for the phases using a solution interval of ~10 s. Following the application of these corrections, a second iteration solved for both amplitude and phase. The flux densities determined for each of the complex gain calibrators observed in the present study are summarized in Table A.1.

Imaging was performed using the AIPS task imagr with robust +1 weighting of the visibilities. Table A.2 provides a summary of the resulting properties of the image and the synthesized beam. Flux densities for each of the target stars were then determined based on Gaussian fits to these images using the AIPS task jmfit). All of the stars were spatially unresolved, hence the measured peak flux density is taken to represent the total flux.

3.4. Visual spectroscopy

Visual domain spectra for the program stars were obtained from the following archives:

• unpublished Heros data, see Štefl &Rivinius (2005);

• the ELODIE archive described by Moultaka et al. (2004);

• the ESPaDOnS archive at CADC⁷. ESPaDOnS is mounted at the CFHT on Hawaii and described by Donati et al. (2006);

• the TBL/Narval archive⁸. NARVAL is a copy of the ESPaDOnS instrument, adapted to the 2 m TBL;

• the BeSS database (Neiner et al. 2011). However, data from this source were not used for velocity measurements; only for equivalent width (EW) measurements and general trend descriptions.

3.5. Visual polarimetry

Visual spectropolarimetric measurements were taken from the archives of the HPOL spectropolarimeter⁹ mounted on the Pine Bluff Observatory (PBO) 36” telescope (Meade et al. 2012).

4. SED modeling

4.1. Model description

The adopted model of a Be star system consists of a combination of a fast rotating central star and a VDD.

The central star is approximated by a spheroidal, rotationally oblate shape that is gravity darkened according to the von Zeipel law (von Zeipel 1924). Its surface is divided into a number of latitude bins, and each one is assigned with a model spectrum (Kurucz 1979) corresponding to the local values of T_eff and log g.

The rotating central star is thus completely described by four independent parameters: mass M, polar radius R_p, luminosity L, and rotation rate W. The W parameter is defined as W = v_rot/v_orb, where v_rot is the rotational velocity at the stellar surface and v_orb is the Keplerian orbital velocity directly above the stellar surface. The gravity darkening exponent β that appears in the classical form of the gravity darkening law, T_eff ∝ g^β (von Zeipel 1924), is a function of W as given by Espinosa Lara & Rieutord (2011).

The VDD model is adopted in a simplified parametric form, in which the density structure follows a power law decrease in the radial direction and has a hydrostatic equilibrium structure in the vertical direction: $\mathit{\rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{e}}}\right)}^{\mathrm{-}\mathit{n}}\exp \left(\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{H}}^{\mathrm{2}}}\right)\mathit{,}$(2)where r and z are respectively the radial and vertical cylindrical coordinates, ρ₀ is the density of the disk close to the stellar surface, R_e is the equatorial radius of the rotationally deformed central star, n is the density exponent, and H is the disk scale height. The disk flares with the following dependence of the disk scale height H on the distance r from the star: $\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{,}$(3)where ${\mathit{H}}_{\mathrm{0}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}{\mathit{v}}_{\mathrm{orb}}^{-1}{\mathit{R}}_{\mathrm{e}}$ is the pressure-supported scale height at the base of the disk, c_s = [ (k_BT_k) / (μm_a) ] ^{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_a is the atomic mass unit. T_k is the parameter that determines the scale height, and we choose it to correspond to 60% of the effective temperature at the stellar equator, following Carciofi & Bjorkman (2006).

For n = 3.5, this structure represents the subsonic parts of an isolated, isothermal VDD in a steady-state (when mass from the central star was fed to the disk at a constant rate for a sufficiently long time) and composed purely of hydrogen (Bjorkman & Carciofi 2005). However, there are effects that can alter the value of n. Local non-isothermalities in the disk can affect both the flaring exponent and the radial density fall-off (Carciofi & Bjorkman 2008), while a non-steady mass feeding rate (Haubois et al. 2012), and the accumulation effect caused by a close orbiting companion (Panoglou et al. 2016) can also affect the density profile in quite complicated ways. For instance, Haubois et al. (2012) found that a forming disk (i.e., a disk being actively fed by the star but that has not yet reached a steady-state configuration) has a sharp density profile, which, if approximated by a power-law, results in n larger than 3.5; conversely, a dissipating disk (i.e., a disk that is no longer being fed and is therefore accreting back to the central star) has a much flatter density exponent (n ≈ 3.0). Density exponents smaller than 3.5 are also believed to ensue from the accumulation effect in Be binaries. In this regard, varying n in our models allows us to accommodate these processes in the studied targets.

The only remaining free parameter of the disk is its outer radius R_out. If not set to a very high number corresponding to a steadily fed isolated disk, this parameter roughly corresponds to the truncation radius caused by the tidal influence of a possible binary companion. In this work, we adopt a sharp truncation, in which no material exists past R_out. Detailed hydrodynamic calculations, however, show that this is not a realistic assumption, because some material spills over past the truncation radius. This material may still contribute to the disk observables (Panoglou et al. 2016).

To compare the model with observations, we additionally need the value of interstellar reddening E(B − V) and two geometrical parameters: the inclination angle i (0° for pole-on orientations) and the distance d.

4.2. Modeling procedure

To calculate the synthetic observables, we employ the code HDUST (Carciofi & Bjorkman 2006, 2008). All the objects are modeled as a combination of a fast-rotating central star surrounded by a VDD. The code uses the Monte Carlo method to solve the radiative transfer, radiative equilibrium, and statistical equilibrium in 3D geometry for the given density and velocity distributions for pure hydrogen gas. For details on the hydrodynamics of VDDs, we refer to Carciofi (2011) and Krtička et al. (20111).

The modeling procedure is similar to the one tested by Klement et al. (2015) when using a multi-technique and multi-wavelength dataset. In the present work we used only the SED flux measurements to constrain the models. As this limits the constraining ability of the modeling, we keep several parameters fixed, specifically the stellar mass, rotation rate, inclination and distance, and we focus instead on determining the parameters that are best constrained by the SED structure, namely the interstellar reddening, the stellar radius and luminosity, the disk parameters (ρ₀, n) and the outer disk radius.

The multi-epoch data were combined into a single dataset, as the variability of the studied stars was found to be reasonably low. Where signs of variability are present, the combined data represent average properties of the disks in the last decades (for details see the results for individual targets in Sect. 5).

We started by searching the literature for spectral type classifications of the programme stars. The spectral types were used to fix the masses of the central stars as given by Harmanec (1988) for the fundamental parameters of main sequence B-type binaries. The rotation rate was kept at the typical value of W = 0.8 (Rivinius et al. 2013). This value corresponds to a R_e/R_p ratio of 1.32. The gravity darkening exponent was correspondingly set to 0.1728, as given by Espinosa Lara & Rieutord (2011). The distances to the stars were taken from the Hipparcos measurements (van Leeuwen 2007), and the inclinations from interferometric studies, when available. For the shell stars, we assumed i > 70°. The adopted parameters are listed in Table 3.

Interstellar reddenning presents an additional model parameter E(B − V), which needs to be determined in order to be able to compare the observations to the model. We used the 2200 Å absorption bump (caused by the interstellar medium) to determine the E(B − V) of the programme stars (we refer to Zagury 2013, for discussion of the origin of the bump). During the procedure the observed UV spectra are dereddened using the reddening curve of Fitzpatrick (1999) with E(B − V) as a free parameter, until the 2200 Å bump disappears and the spectrum becomes flat in that region. For the extinction R_V we used the standard value of 3.1 (Savage & Mathis 1979). Although values other than 3.1 may be observed in nebular and star-forming regions, for nearby field stars such deviations are less likely. The uncertainties of the determined values of the E(B − V) parameter were inferred by a Markov Chain Monte Carlo (MCMC) method provided by the emcee Python module¹⁰ (Foreman-Mackey et al. 2013).

We proceed to the remaining physical parameters of the central star, namely the polar radius R_p and luminosity L. As discussed in Klement et al. (2015), for non-shell stars or shell stars with tenuous disks, the disk has very little influence on the observed UV spectrum. In such cases and within a reasonable parameter range, changing R_p (which in effect means changing the effective temperature as well) influences the slope of the model UV spectrum only, while changing L influences mainly its level. This allows us to constrain these two parameters independently by computing a grid of purely photospheric, fast-rotating models (similar to Klement et al. 2015). Furthermore, for tenuous disks, the visual and near-IR parts of the SED are also devoid of disk influence. The wavelength interval used to determine R_p and L was therefore selected on a case-by-case basis (we refer to Sect. 5 for details on individual targets).

The best-fit R_p and L parameter combination was determined via χ² minimization, where for the best-fit model ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}{\mathit{\chi }}_{\min }^{\mathrm{2}}$. Since we were dealing with SED modeling, where both the wavelengths and fluxes span many orders of magnitude, we used the logarithmic version of χ²: ${\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\left(\frac{\log {\mathit{F}}_{\mathrm{obs,}\mathit{i}}\mathrm{-}\log {\mathit{F}}_{\mathrm{mod,}\mathit{i}}}{{\mathit{F}}_{\mathrm{err,}\mathit{i}}\mathit{/}{\mathit{F}}_{\mathrm{obs,}\mathit{i}}}\right)}^{\mathrm{2}}\mathit{,}$(4)where F_obs is the observed flux, F_mod is the model flux, F_err is the observed flux error, and N is the total number of fitted data points. The uncertainties of the best-fit parameters were estimated from the $\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{0.683}}^{\mathrm{2}}\mathrm{=}{\mathit{\chi }}_{\mathrm{mod}}^{\mathrm{2}}\mathrm{-}{\mathit{\chi }}_{\min }^{\mathrm{2}}$ contour corresponding to the 68.3% confidence level (1σ), where $\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{0.683}}^{\mathrm{2}}$ is a function of the number of degrees of freedom of the fit (we refer to, e.g., Chap. 11 of Bevington & Robinson 1992). In some cases the parameter uncertainties turned out to be lower than the step size adopted for the parameter grid. Since the quality of the model fit was already sufficient for the aims of this work, we did not put in the computational effort to further fine-tune the model parameters and thus lower the uncertainty estimates, and we instead conservatively selected the corresponding step size as the parameter uncertainty.

The disk itself is non-isothermal with a temperature structure computed by HDUST, but its density structure is not solved self-consistently with the temperature solution and remains the same throughout the simulation. For such a simplified VDD model we need only three free parameters to fully describe it: the base density ρ₀, the density exponent n, and the disk outer radius R_out. We set an upper limit to the R_out parameter, so that it cannot exceed the value of the critical radius R_c estimated by means of Eq. (1), in which the sound speed c_s was computed using the mass-averaged temperature of the disk.

In our model, we assume a single value of n throughout the whole disk and by doing so we are neglecting local non-isothermal density effects in the disk. VDD theory predicts a temperature minimum to occur in the inner parts of the disk due to its high opacity to the stellar radiation. However, since we are interested in the overall density structure of the VDD, of which the inner non-isothermal part is only a small fraction, these effects are not important for our analysis.

The modeling procedure was executed as follows. First, we fixed R_out to an intermediate value of 50 stellar radii, and a grid of models with ρ₀ and n in the range 1× 10^-13–1× 10^-10 g cm^-3 and 1.5–5.0, respectively, was computed. To determine the parameters ρ₀ and n, we fit the model grid to the part of the SED where the disk contribution is significant, but not to the long wavelength radio data, that are sensitive to the disk size. While dense disks start to contribute to the SED already at visual wavelengths, excess radiation due to a tenuous disk might become observable only at mid-IR wavelengths. Similarly, the truncation effects in very dense disks will influence the SED at shorter wavelengths as compared to less dense disks. Therefore, the SED region used to constrain ρ₀ and n will be different from star to star and will shift towards shorter wavelengths with increasing density of the disk (for details on individual targets see Sect. 5). The uncertainties of the best-fit values of the parameters ρ₀ and n were determined in the same way as for the central star parameters R_p and L, described above.

The modeling procedure gets more complicated for shell stars with dense disks. In these cases, the disk reduces the UV flux (and sometimes also portions of the visual and near-IR flux) emitted by the central star and the observed UV spectrum cannot, therefore, be used to constrain R_p and L in a straightforward way. We adopted an iterative modeling procedure for such stars: A first look estimate for R_p and L was done using the photospheric models only, after which the disk was included in the model to get the first estimates for the disk parameters ρ₀ and n from the IR fluxes. With the values obtained in this way, the magnitude of the influence of the disk on the UV spectrum could be estimated and taken into account for the next iteration of photospheric parameters. After a few such iterations, a satisfactory fit to both the UV and IR parts of the SED was obtained.

After determining the best-fit values of ρ₀ and n, we can proceed to the determination of the R_out parameter. As mentioned above, the truncation effects will become detectable at progressively shorter wavelengths with increasing density of the disk. However, with disk sizes ≳20 R_e, the influence does not fall shortwards of sub-mm/mm wavelengths. Therefore, the mm and cm fluxes are suitable for constraining the sizes of dense disks, while cm fluxes only constrain the disk size of the tenuous ones (e.g., β CMi). In the latter case, the sub-mm/mm region was included when searching for the best-fit values of ρ₀ and n. Measurements corresponding to upper limits were used for visual checking, but not for the modeling procedure itself. The uncertainty of the R_out parameter was obtained in the same way as for the previously mentioned model parameters.

5. Results

In this section we first review the results for β CMi from Klement et al. (2015). Then we present the results for stars for which we have new JVLA measurements (η Tau, EW Lac, ψ Per, γ Cas), followed by the results for β Mon A, which has new sub-mm data from APEX/LABOCA. We also give a brief overview of what can be learned from the available visual spectroscopy and polarimetry for each target. The derived model parameters for all stars are listed in Table 4. In Table 4 we also include the values of R_c calculated using Eq. (1).

Fig. 1 β CMi. The best-fit model is plotted as the black solid line, the photospheric contribution as the dashed black line, and a non-truncated disk model (Rout = Rc) as the black dotted line. Upper: SED in the 1000–10 000 Å interval. Middle: SED from 0.1 μm to 10 cm. The inset shows model predictions for different disk sizes. The observed fluxes are plotted with error bars, which are smaller than the symbol size in most cases. Lower: residuals structure of the best-fit model with Rout = 40 Re. The photospheric and non-truncated models are also plotted in relation to the best-fit model, highlighting where these models deviate from the best-fit one.

5.1. β CMi (HD 58715; HR 2845)

The star β CMi was studied in extensive detail by Klement et al. (2015). As in the present study, the HDUST code was used to model the SED from the UV to the radio, but additionally a large data set of spectroscopic, polarimetric, and interferometric observations was used to further constrain the model. Here we present a slightly adjusted model of the SED produced using an updated version of the HDUST code, which fixed some numerical problems that were slightly affecting the correct computation of synthetic observables. More details on the fixes of the code and the changes regarding the results of modeling of the observables besides the SED are given in Klement et al. (2017).

We recomputed the model and allowed the ρ₀, n and R_out parameters to vary in a small interval around the values originally determined by Klement et al. (2015), while keeping the central star parameters the same as in the previous study. The best-fit value of base density ρ₀ remains the same, but a significant change to the model concerns the density exponent n: the SED up to mm wavelengths is now reproduced exceptionally well when using a value of n = 2.9. We recall that in the original study, the very inner parts of the disk were compatible with a steeper (n = 3.5) density fall-off, while in the remainder of the disk it followed a shallower profile (n = 3.0).

Regarding the outer parts of the disk, the original conclusion was that the disk is truncated at a distance of $\mathrm{3}{\mathrm{5}}_{-5}^{\mathrm{+}\mathrm{10}}$ R_e, as revealed by the observed radio flux at 2 cm. Here we include three additional VLA measurements at cm wavelengths that are available from the literature. The revised model and additional cm data led to a revised best-fit disk size of $\mathrm{4}{\mathrm{0}}_{-5}^{\mathrm{+}\mathrm{10}}$ R_e. The best-fit model is plotted in Fig. 1 as a solid black line, together with the purely photospheric flux (dashed black line) and a non-truncated disk model (dotted black line). The model reproduces the full SED almost perfectly, with the exception of the data point at 6 cm, which lies slightly above our best-fit model.

5.2. η Tau (Alcyone; HD 23630; HR 1165)

This star from the Pleiades cluster has been mostly classified in the literature as a B7IIIe star (Slettebak 1982; Lesh 1968; Hoffleit & Jaschek 1991). The inclination is apparently of an intermediate value. In previous interferometric studies it was quoted to be higher than 18° (Quirrenbach et al. 1997) and equal to 41° (Tycner et al. 2005). We adopt the latter value for our model.

No companion was revealed in adaptive optics observations (Roberts et al. 2007). A total of 47 spectra suitable for RV measurement were found in the archives. The RV was measured in the Mg ii 4481 line. The mean RV is 6.9 ± 2.1 km s^-1. Heros measurements around MJD 52000 are on average approximately 4 km s^-1 lower than ESPaDOnS measurements at MJD 55 000, but given the different instrument parameters, the values are well compatible and within 3σ of one another. However, we note that this by no means excludes the possible presence of a binary companion. The IR Ca triplet emission, indicative of binarity, is not present in the spectra.

The spectroscopic data show that the Hα emission was at a constant level of I/I_c = 2.2, and W_λ = − 3.5 Å until approximately 2005, and since then steadily decreased by a small but detectable amount to I/I_c = 2.0, and W_λ = − 2.5 Å. The spectropolarimetric measurements available from the HPOL database show no change of polarization degree or angle between 1992 and 1999.

From the comparison of IRAS and AKARI/WISE measurements, which are separated in time by around two decades, we also see that the variability of the disk has been relatively low. This is a common feature of late-type Be stars. The comparison of the observed mm fluxes from JCMT and IRAM shows signs of slight variability at mm wavelengths on the time scale of years. The archival VLA and the new JVLA measurements are very much consistent with one another, with the exception of the observed fluxes at 3.6 cm, which however still lie within 2σ of one another.

Fig. 2 As in Fig. 1, but for η Tau.

Fig. 3 As in Fig. 1, but for EW Lac.

Fig. 4 As in Fig. 1, but for ψ Per.

The results of the SED modeling of η Tau are plotted in Fig. 2. The non-shell nature of η Tau allows us to use the UV spectrum to determine the central star parameters R_p and L without any complications. Moreover, the disk of η Tau is apparently very tenuous and starts to contribute to the observed SED only at around 5 μm. We thus use both the UV spectrum and the visual photometry to constrain R_p and L. Changing the disk size influences the cm fluxes only and has a negligible influence on the mm fluxes (as expected for tenuous disks). Therefore, the SED interval from 5 μm to mm was used to find the best-fit values of the disk parameters ρ₀ and n.

The observed radio fluxes are lower than expected from a non-truncated VDD model. The disk size best reproducing the cm data is ${\mathit{R}}_{\mathrm{out}}\mathrm{=}\mathrm{4}{\mathrm{0}}_{-5}^{\mathrm{+}\mathrm{10}}$ R_e. We interpret this as a clear sign of truncation of the outer disk. However, a simply truncated VDD does not offer a good fit either, as the observed SED slope in the radio is flatter than the model predicts. This is most easily seen in the residuals of the best-fit model with respect to the radio measurements (bottom panel of Fig. 2).

5.3. EW Lac (HD 21750; HR 8731)

EW Lac was classified as a B3IVe-shell star (Slettebak 1982; Rivinius et al. 2006), and a B4IIIep star (Lesh 1968; Hoffleit & Jaschek 1991). No previous interferometric determination of the disk inclination exists in the literature. However, the shell nature of the star restricts the inclination angles to ≥70°. For our study we adopt the value of i = 80°, which was found to reproduce well the shapes of the Balmer and Paschen jumps.

EW Lac has shown remarkable variations in its shell line appearance in the past, most notably from 1976 to 1986, and a similar episode lasted from 2007 until 2012 (Mon et al. 2013). The variations are best explained in terms of V/R variability, that is, changes in the peak height ratio of the violet and red side of the emission. Small V/R differences were present in almost all spectra at our disposal, taken from 1998 to 2014. However, beginning in the 2007 to 2014 period, that is, for the last V/R cycle, the amplitudes were greatly enhanced.

There are 13 suitable spectra to measure RVs, from which the presence of a companion could be neither confirmed, nor excluded. Values measured in the Fe ii 5169 line obtained from Heros spectra, from MJD 51 000 to 52 500, are at approximately − 22 km s^-1, while ESPaDOnS values, obtained around MJD 55 400, are at − 28 km s^-1. Although this change is very likely real, it is probably due to a change in the global density oscillation of the disk rather than binarity. Also, changes in the Hα equivalent width are due to the variable strength of the central absorption core, rather than changes in the strength of the emitting peaks.

Like for η Tau, the HPOL data, taken from 1989 to 1995, do not show any change in polarization degree or position angle; neither is the infrared Ca triplet observed in emission in our spectra.

Comparison of the IRAS and AKARI/WISE measurements reveals a variation in the IR excess emission that is somewhat stronger than in the other targets in our sample. The mm fluxes observed by JCMT and IRAM are consistent within the respective error bars. The comparison of the VLA and JVLA measurements also reveals a slight discrepancy in that one of the VLA measurements lies a few σ above the JVLA fluxes. The two VLA observations at 2 cm are slightly inconsistent, as the upper limit measurement lies below the detection. This suggests a possible mild variability of cm fluxes on the time scale of a year.

The results of the SED modeling of EW Lac are plotted in Fig. 3. Due to the shell nature of the star and the high value of the disk density, the observed UV spectrum of EW Lac is slightly dimmed by the surrounding disk. The disk excess starts to be significant at a wavelength of around 1 μm, allowing for the inclusion of the visual part of the SED for determining the central star parameters. Since the disk is dense, the IR part of the SED was used to determine ρ_O and n, while the mm and radio part, sensitive to varying the disk size, was used to determine R_out. The resulting SED fit up to far-IR is satisfactory.

The radio observations, including the mm measurements, are inconsistent with a non-truncated disk. The resulting disk size for the case of simple truncation is $\mathrm{10}{\mathrm{0}}_{-20}^{\mathrm{+}\mathrm{30}}$ R_e. Moreover, the new JVLA data set reveals a similar result as given above for η Tau – the observed slope of the SED in the radio is not as steep as would be expected for a simply truncated disk.

5.4. ψ Per (HD 22192; HR 1087)

The star ψ Per has been classified as a B5Ve-shell star (Lesh 1968; Hoffleit & Jaschek 1991), a B4IVe star (Underhill 1979) and also as a B5IIIe star (Slettebak 1982; Rivinius et al. 2006). The inclination angle is >62° according to Quirrenbach et al. (1997) and 75°±8° according to Delaa et al. (2011); for this study we adopt the latter value.

No binary companion is evident (Mason et al. 1997; Delaa et al. 2011). Yet, for this star, the infrared Ca triplet is observed in clear emission in all available spectra covering that region. Like for EW Lac, the shell nature permits the measurement of the RV with a precision of a few dozen m s^-1. While the RV is varying by approximately 1500 m s^-1 peak to peak, there are not enough spectra to say whether this is due to a companion, that is, clearly periodic, or due to small changes in the disk, in which case changes could be merely cyclic or secular.

For ψ Per, the HPOL data show a clear increase of polarization degree from 0.4% to 0.8% from 1992 to 1994, which then remains at 0.8% until at least 2000. All available spectra were taken in 2000 to 2014, that is, after the polarization changed, and show a strong and stable disk, with I/I_c = 6 and W_λ = − 41 Å. The V/R ratio of the emission lines is unity; apparent deviations observed, for example, in Fe ii 5169, are due to blends.

No variability in the IR excess is revealed by comparing the IRAS and AKARI/WISE data. The mm measurements from JCMT and IRAM are mostly consistent, although one of the JCMT detections lies slightly above the others. The old VLA and the new JVLA data are in good agreement, although the upper limit of Apparao et al. (1990) at 2 cm lies below the subsequent detections. It is unclear whether this is caused by radio variability or by the low signal-to-noise ratio of the observations.

The results of the SED modeling of ψ Per are plotted in Fig. 4. Although ψ Per is a shell star, the disk base density is only of an intermediate value and does not significantly influence the UV and visual part of the SED. Therefore the SED interval up to 1 micron was used to constrain the central star R_p and L. The disk contribution to the observed SED starts to be significant at near-IR wavelengths. Although the value of ρ₀ is only intermediate, the fact that n is very low makes the disk dense in its outer parts. Correspondingly, the model SED was found to be sensitive to changing the disk size already at mm wavelengths, therefore it was only the IR portion of the spectrum that was used to constrain the disk density parameters. The resulting fit from UV to far-IR wavelengths is exceptionally good.

The observed radio SED slope is once again inconsistent with both a non-truncated and a simply truncated disk. The resulting disk size best reproducing the whole radio data set is $\mathrm{10}{\mathrm{0}}_{-15}^{\mathrm{+}\mathrm{5}}$ R_e. The result again suggests that the disk density fall-off regime changes to a steeper one in the outer parts, but the disk is not simply cut-off, as is assumed in our model.

To this date, ψ Per is the only star that has been angularly resolved in radio (λ = 2 cm) along the semi-major axis of its disk (Dougherty & Taylor 1992). In that study, the Gaussian fit to the azimuthally averaged visibility data indicated a best-fit disk size (FWHM) of 74 mas, corresponding to 392 R_e. In Fig. 5 we plot the original visibility curve of Dougherty & Taylor (1992), along with the curves derived from our models with different disk sizes (azimuthally averaged). Since the absolute flux density values for each of our models differ, we scale the curves so that at zero baseline the visibility is equal to unity. The visibility curve of our best-fit model (R_out = 100 R_e) shows that such a disk would not have been resolved by the observations of Dougherty & Taylor (1992), although the model flux density at 2 cm (1190 μJy) is relatively close to the measured value (804 ± 60 μJy). The model with the disk size of 500 R_e shows almost perfect agreement with the shape of the observed visibility curve. However, the flux density of this model (9.31 mJy) overestimates the measured value by more than a factor of 10. As for the apparent disk sizes of our models, 2D Gaussian fits give FWHMs of 29.5 × 20.1 mas and 85.2 × 47.8 mas for the best-fit and the 500 R_e disk models, respectively. The disagreement between our model and the angular size and fluxes of Dougherty & Taylor (1992) cannot be solved by the simple models that we employ, and will be the subject of future research.

Fig. 5 ψ Per: azimuthally averaged visibility data of Dougherty & Taylor (1992, black points with error bars) and a Gaussian fit to them (thick black line, FWHM of 74 mas corresponding to 518 Re) over-plotted with the visibility curves derived from our azimuthally averaged models with several different disk sizes indicated in the legend. The dotted line shows the visibility curve of a point source (visibility equal to unity with increasing baseline length).

5.5. γ Cas (HD 5394; HR 264)

The star γ Cas was one of the first two Be stars to be discovered (Secchi 1866). It was identified as an X-ray source and there is an ongoing debate as to the origin of the puzzling X-ray emission properties (e.g., Motch et al. 2015). γ Cas is a known single-line spectroscopic binary. The set of orbital parameters was recently revised by Nemravová et al. (2012), who concluded that the orbit is circular with an RV semi-amplitude of 4.30 ± 0.09 km s^-1 and a period of 203.52 ± 0.08 days. γ Cas has a very complicated observational history, as it has appeared as a shell star, briefly as a normal B star and since a few decades ago as a Be star seen under intermediate inclination (see Harmanec 2002, for details).

γ Cas has been classified in the literature as a B0.5IVe star (Slettebak 1982). For the inclination, we adopt the value of 42°, found in a recent interferometric study by Stee et al. (2012).

In agreement with the suspicion that the infrared Ca triplet may be linked to binarity (Koubský et al. 2011), γ Cas does indeed show emission in these spectral lines. The strength is somewhat lower than that of the Paschen lines that the Ca triplet lines are blended with, but clearly discernible. The shape of the Hα line also seems to indicate binarity, since it often shows no clear peaks or symmetry, which are signs of disturbances in the disk.

Fig. 6 As in Fig. 1, but for γ Cas.

Judging from the Hβ profiles, which behave closer to the expectations for a classical Be star than the Hα profiles, γ Cas underwent V/R variability before 2000, which weakened between 2000 and 2005, and is since absent. The available HPOL data, spanning 1991 to 2000, show a constant polarization. Since γ Cas is the brightest Be star in the Northern hemisphere, there is a great wealth of spectroscopic data available from the archives. Here spectra taken between 1996 and 2014 are considered. The value for I/I_c of Hα varied, but remained in the range between approximately 4 and 5. In other words, while the disk of γ Cas is not stable in the same sense as the one of, for instance, ψ Per, it is clearly not currently in a stage of secular build-up or decay, but fluctuating around a certain state.

The IR measurements do however reveal signs of variability, as the IRAS fluxes are several σ below those of AKARI/WISE. The same is true for the observed mm fluxes, as the JCMT and IRAM measurements are discrepant by several σ. This suggests mild variability in the disk, consistent with what is observed in the Hα emission. The VLA and JVLA measurements are in very good agreement, except for the 3.6 cm VLA upper limit which lies slightly below the JVLA detection.

The results of the SED modeling of γ Cas are plotted in Fig. 6. The disk contribution to the observed SED is apparent already at visual wavelengths, therefore we use only the UV spectrum to constrain the central star parameters R_p and L. Owing to the high density of the disk, the observed fluxes from optical to far-IR were used to determine the disk parameters ρ₀ and n. The overall model fit to the observed SED from UV to far-IR region is satisfactory with the exception of the visual and near-IR fluxes, which are overestimated by the model.

The result from the modeling of the radio SED is that the disk is clearly truncated, similar to the previous targets. The disk size best reproducing the radio data is $\mathrm{3}{\mathrm{5}}_{-5}^{\mathrm{+}\mathrm{5}}$ R_e. However, the observed slope of the radio SED is again flatter than the slope of the truncated model.

Fig. 7 As in Fig. 1, but for β Mon A.

5.6. β Mon A (HD 45725; HR 2356)

β Mon A is a component of a visual multiple system, with the B and C components separated by 7.1 and 10 arcsec, respectively (Taylor et al. 1990). At such distances the companions cannot have any tidal influence on the disk of the A component (Rivinius et al. 2006). The star has been classified as a B3Ve star (Lesh 1968; Hoffleit & Jaschek 1991), B4Ve-shell star (Slettebak 1982; Rivinius et al. 2006), and a B2III star (Maranon di Leo et al. 1994). Strong V/R variations with a period of 12.5 years have been observed in the past (1930–1960, Taylor et al. 1990).

The shell nature of the star restricts the inclination angle to be close to edge-on. Unfortunately no interferometric studies of this star have been published (although interferometric measurements do exist and are public - Rivinius et al. 2015); the only reference for the inclination is 67° (Frémat et al. 2005). We adopt an inclination of 70° for our model.

There are 48 available spectra taken in four observing seasons between MJD 51 000 and 56 000. The narrow circumstellar line Fe ii 5169 shows small RV variations with an amplitude of approximately 3 km s^-1. Again, however, it is not clear whether this variation is due to binarity or changes in the disk, since β Mon A is a V/R variable Be star.

No HPOL data are available for β Mon A. Judging from the line emission, the disk is generally strong, with values of I/I_c between 5 and 6. It also underwent a weak V/R cyclic variability in these years. The spectra do show IR Ca triplet emission on a similar level to that of γ Cas.

The disk properties as shown by the observed SED seem to have been relatively stable in the last decades, with only the AKARI measurements being slightly discrepant with respect to those from IRAS and WISE. The sub-mm/mm measurements indicate a slight variability, similar to the previous targets. New JVLA observations are not available for this star and the data set at cm wavelengths consists only of a single VLA detection at 2 cm and three upper limit measurements at 2, 3.6, and 6 cm. The upper limit at 2 cm lies a few σ below the detection, suggesting a mild variability of the radio fluxes.

The results of the SED modeling of β Mon A are plotted in Fig. 7. β Mon A is a shell star and the disk influence is already clear in the UV spectrum, which is slightly dimmed in comparison with a purely photospheric spectrum. The disk of β Mon A is dense and significant excess radiation is present already at visual wavelengths. Only the IR part of the SED was used to search for the disk parameters, as the disk has a high density. The resulting model SED agrees with the observations reasonably well.

For β Mon A the signs of the disk truncation are the weakest among our targets. However, the APEX/LABOCA and JCMT fluxes lie more than 3σ below what is expected of a non-truncated disk. The upper limit measurements at 2 and 3.6 cm are also inconsistent with a non-truncated disk. Nevertheless, additional radio measurements are necessary to better constrain the disk size and confirm the signs of truncation. However, even with the available data, the best-fit disk size is 110 ± 40 R_e, which is significantly lower that the critical radius value of 260 R_e.

6. Discussion

The VDD model in its parametric form reproduces the observed SED of the studied stars generally well. The resulting values of the density exponent n are in all cases lower than the canonical value of 3.5 derived for flaring isothermal VDDs in steady-state. A possible explanation for this was recently given by Vieira et al. (2017), who studied the IR SEDs of 80 Be stars and determined their ρ₀ and n parameters using a semi-analytic model based on the pseudo-photosphere concept (Vieira et al. 2015). The determined values of n are generally not equal to 3.5, with the majority of the studied disks showing n between 2.0 and 3.0. These results were proposed to be explained by the dynamical state that we see the disks in. According to the dynamical models (Haubois et al. 2012) and assuming additional effects such as cooling by heavier elements (not taken into account in our model), the disks with n ≥ 3.5 should correspond to disks that are in the process of formation, disks with n between approximately 3.0 and 3.5 are in a steady-state, and disks with n ≲ 3.0 are in the process of dissipation. As presented in Fig. 7 of Vieira et al. (2017), approximately two thirds of the observed Be disks were found to be in the dissipating phase, suggesting that the time scales for disk dissipation are longer than those for disk formation.

Over-plotting the results of the present study on Fig. 7 of Vieira et al. (2017) shows that the log ρ₀ and n parameter combinations determined for the stars in our sample fall within the two central contours corresponding to the probability density of 0.06 (Fig. 8). However, contrary to the study of Vieira et al. (2017), the majority of our targets seem to be in steady-state (4 out of 6). This apparent disagreement (which might not be significant given the small number of stars in our sample) is probably due to the fact that Waters et al. (1991) favored Be stars with large infrared and radio excesses in their study. Another apparent inconsistency is that the two targets for which n is well below 3.0 (EW Lac and ψ Per) and should, according to Vieira et al. (2017), be in a dissipation state, also seem to have stable disks over long time-scales. We speculate that if these disks are truncated by unseen binary companions (see below), the shallower density profile might be due to the accumulation effect, which causes the disk density profile to become shallower than the steady-state value (Panoglou et al. 2016).

An SED turndown (i.e., a steepening of the SED slope) is clearly observed between the far-IR and radio wavelengths in five out of the six studied stars. The one exception is β Mon A, for which the presence of an SED turndown is inconclusive.

The results of the modeling show that, to a first approximation, the observed SED turndown can be explained by assuming truncated disks. We recall that there are two physical effects that are predicted to affect the outer disk structure, that could be the cause of the SED turndown: the tidal influence from a close binary companion, or the transonic transition, that occurs when the outflow velocity of the gas in the disk becomes supersonic.

Fig. 8 Results for the ρ0 and n distributions of the sample of Be stars studied by Vieira et al. (2017). The main panel shows the probability density (gray scale contours), while the upper and right panels show the distributions for the individual parameters. The contour values correspond to the probability density levels. Superimposed as filled circles are the results of Vieira et al. for individual stars (colors indicate the stellar effective temperatures computed by Frémat et al. 2005). Over-plotted as yellow star symbols are the results for the stars studied in this work.

We found that the disk truncation occurs at much smaller distances from the central star than where the supersonic point (the critical radius, R_c, in Table 4) is expected to lie. Therefore we conclude that the supersonic regime is an unlikely cause for the SED turndown observed in our sample. As mentioned in Sect. 2, the effects of the transonic transition in a typical Be disk are likely to be observable only at even longer wavelengths. The one exception is EW Lac, for which the non-truncated disk model SED structure in the radio (dotted line in the lower panel of Fig. 3) suggests that for this star (and Be stars with similarly dense disks), the effect of the transonic transition may already be observable at cm wavelengths. However, we note that the formula for the critical radius (Eq. (1)) is only approximate and the resulting values of R_c are therefore uncertain.

The disk sizes determined by our modeling represent sizes of truncated disks that best reproduce the data, but it is clear from the radio SED slope that the situation is not that simple: some matter seems to overflow beyond the truncation radius (which is generally not equal to our parameter R_out), from where it non-negligibly contributes to the radio SED. Investigating the role of the material past the truncation radius on the SED will be the subject of a future study.

For the previously studied star β CMi, the presence of the companion was not confirmed at the time of the analysis of Klement et al. (2015). The disk truncation was nevertheless suggested as being caused by the tidal influence of an unseen binary companion, with further observable evidence, such as the Ca triplet in emission, supporting this option (we refer to Sect. 5 of Klement et al. 2015). A follow-up RV analysis of the Hα line led to the detection of the binary, with an orbital period of 170 days and RV amplitude of 2.25 km s^-1 (Dulaney et al. 2017). The value of the detected period, if a circular, coplanar orbit and a 1 M_⊙ companion are assumed, indicates a semi-major axis of ~50 R_e. Based on the results of Panoglou et al. (2016), the truncation radius is expected to lie near the 3:1 resonance with the binary orbit, corresponding to ~25 R_e. The result that the disk is truncated at ~40 R_e rather favors the 3:2 resonance with the orbit as the location of the truncation radius for the case of β CMi.

For the other confirmed binary in our sample, γ Cas, the semi-major axis of the orbit is a ≃ 1.64 au (Nemravová et al. 2012), corresponding to 38 R_e of our model. However, here the comparison with our best-fit disk size is difficult, since the slope of the radio SED is not reproduced well. This means that if, for example, the cm data were ignored, the mm measurements would point to a much smaller disk size than was derived from the whole radio data set. This is also the case for the remaining objects with new JVLA data (η Tau, EW Lac, and ψ Per).

In conclusion, modeling of the radio SED offers the possibility to indirectly detect previously unknown companions of Be stars. A systematic search for RV periods in stars which show truncation effects (especially shell stars), but for which no companion was reported, is necessary to confirm that it is indeed binary companions that cause the SED turndown by truncating the disk, and not any other so-far undetected mechanism operating in Be star disks. We performed a search for RV variations in the available spectra of the unconfirmed binaries in the sample. However, the data set was limited and the RV amplitudes caused by low-mass companions are likely of the same order as the value 2.25 km s^-1 derived for β CMi for late-type Be stars, and even smaller for earlier types (since the mass ratio is higher). Detection of such RV amplitudes would probably require dedicated high-resolution spectroscopic campaigns and careful analysis similar to the one performed by Dulaney et al. (2017). If the companions are confirmed and are found to be sdO/sdB stars, this may establish a common evolutionary path for Be stars and help confirm a possible cause for Be stars being fast rotators.

We note that the original sample of Waters et al. (1991) was biased towards stars with the flattest spectra, that is, with the least pronounced SED turndowns. The reason is that only stars detected at cm wavelengths were analyzed, while those not detected likely have even more pronounced SED turndowns and therefore stronger signs of truncation.

7. Conclusions

The main conclusions of this work may be summarized as follows.

• The predictions of the VDD model agree very well with theobserved SEDs covering the wavelength interval from the UV tothe far-IR. This further establishes the VDD model as the bestmodel so far to explain the structure of Be disks.

• Comparison of our results with the ones of Vieira et al. (2017) shows that our sample has a much higher prevalence of steeper density profiles (n ~ 3, typically associated with disks in steady-state) than shallower profiles (n ~ 2, usually found in Be stars whose disks are dissipating). However, for these latter cases (namely EW Lac and ψ Per) there is evidence that the disks have remained stable over the last decades, suggesting that the shallow density profile may be caused by the accumulation effect due to the unseen companions (Panoglou et al. 2016) rather than due to the disks being in a dissipation phase.

• A clear steepening of the spectral slope at far-IR to radio wavelengths is observed in the whole sample with the exception of β Mon A, for which the measurements are not conclusive. The effect is present both for confirmed binaries (β CMi, γ Cas) as well as for stars for which no companion has yet been reported.

• We propose that the SED turndown is caused by truncation of the disk. The radio data sets are best reproduced by disks truncated at distances of approximately 30–150 stellar radii from the central star. The most plausible explanation for the disk truncation is the presence of (unseen) binary companions, tidally influencing the outer disk.

• Our simple model for truncation, in which the disk is abruptly cut at the truncation radius, cannot explain in detail the shape of the radio SED. This is in line with hydrodynamical simulations that predict that the disk material overflows past the truncation radius. In a future study we will investigate whether or not this extra material can explain the discrepancies found in this work.

• Comparison of our results with the only previous determination of the disk size in radio (Dougherty & Taylor 1992, ψ Per) revealed unexpected inconsistencies. Our model that best fits the radio flux densities is too small to have been resolved by the VLA. Interestingly, it is the model with the disk size of 500 R_e that shows good agreement with the observed shape of the visibility curve. However, this model clearly overestimates the observed flux density.

Acknowledgments

We acknowledge our recently deceased colleague Stanislav Štefl as being the one who came up with the original idea for this project.

R.K. would like to acknowledge the kind help of Dietrich Baade in the final stages of preparation of the paper.

The research of R.K. was supported by grant project number 1808214 of the Charles University Grant Agency (GA UK).

A.C.C. acknowledges support from CNPq (grant 307594/2015-7) and FAPESP (grant 2015/17967-7).

R.G.V. acknowledges support from FAPESP (grant 2012/20364-4).

D.M.F. acknowledges support from CNPq (grant 200829/2015-7) and FAPESP (grant 2016/16844-1).

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.

This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut fur Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory

The National Radio Astronomy Observatory is operated by Associated Universities, Inc., under cooperative agreement with the National Science Foundation.

This publication makes use of VOSA, developed under the Spanish Virtual Observatory project supported from the Spanish MICINN through grant AyA2011-24052.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

The VLA data presented here were obtained as part of program VLA/10B-143 (AI141). LDM acknowledges support from an award from the National Science Foundation (AST-1516106).

References

Appendix A: Additional tables

All Tables

All Figures

Fig. 1 β CMi. The best-fit model is plotted as the black solid line, the photospheric contribution as the dashed black line, and a non-truncated disk model (Rout = Rc) as the black dotted line. Upper: SED in the 1000–10 000 Å interval. Middle: SED from 0.1 μm to 10 cm. The inset shows model predictions for different disk sizes. The observed fluxes are plotted with error bars, which are smaller than the symbol size in most cases. Lower: residuals structure of the best-fit model with Rout = 40 Re. The photospheric and non-truncated models are also plotted in relation to the best-fit model, highlighting where these models deviate from the best-fit one.

In the text

	Fig. 2 As in Fig. 1, but for η Tau.
In the text

	Fig. 3 As in Fig. 1, but for EW Lac.
In the text

	Fig. 4 As in Fig. 1, but for ψ Per.
In the text

Fig. 5 ψ Per: azimuthally averaged visibility data of Dougherty & Taylor (1992, black points with error bars) and a Gaussian fit to them (thick black line, FWHM of 74 mas corresponding to 518 Re) over-plotted with the visibility curves derived from our azimuthally averaged models with several different disk sizes indicated in the legend. The dotted line shows the visibility curve of a point source (visibility equal to unity with increasing baseline length).

In the text

	Fig. 6 As in Fig. 1, but for γ Cas.
In the text

	Fig. 7 As in Fig. 1, but for β Mon A.
In the text

Fig. 8 Results for the ρ0 and n distributions of the sample of Be stars studied by Vieira et al. (2017). The main panel shows the probability density (gray scale contours), while the upper and right panels show the distributions for the individual parameters. The contour values correspond to the probability density levels. Superimposed as filled circles are the results of Vieira et al. for individual stars (colors indicate the stellar effective temperatures computed by Frémat et al. 2005). Over-plotted as yellow star symbols are the results for the stars studied in this work.

In the text