Revealing the structure of the outer disks of Be stars

TL;DR

This study investigates the outer disk structure of Be stars using radio and multi-wavelength observations, testing the viscous decretion disk model's ability to explain observed spectral energy distribution features, including the radio turndown.

Contribution

It provides the first detailed radiative transfer modeling of the outer regions of Be star disks, incorporating binary effects to explain the radio spectral turndown.

Findings

VDD models reproduce IR to far-IR SEDs well

Radio SED turndown is not fully explained by VDD models

Binary interactions may influence outer disk structure

Abstract

Context. The structure of the inner parts of Be star disks (20 stellar radii) is well explained by the viscous decretion disk (VDD) model, which is able to reproduce the observable properties of most of the objects studied so far. The outer parts, on the ther hand, are not observationally well-explored, as they are observable only at radio wavelengths. A steepening of the spectral slope somewhere between infrared and radio wavelengths was reported for several Be stars that were previously detected in the radio, but a convincing physical explanation for this trend has not yet been provided. Aims. We test the VDD model predictions for the extended parts of a sample of six Be disks that have been observed in the radio to address the question of whether the observed turndown in the spectral energy distribution (SED) can be explained in the framework of the VDD model, including recent theoretical development for truncated Be disks in binary systems. Methods. We combine new multi-wavelength radio observations from the Karl. G. Jansky Very Large Array (JVLA) and Atacama Pathfinder Experiment (APEX) with previously published radio data and archival SED measurements at ultraviolet, visual, and infrared wavelengths. The density structure of the disks, including their outer parts, is constrained by radiative transfer modeling of the observed spectrum using VDD model predictions. In the VDD model we include the presumed effects of possible tidal influence from faint binary companions. Results. For 5 out of 6 studied stars, the observed SED shows strong signs of SED turndown between far-IR and radio wavelengths. A VDD model that extends to large distances closely reproduces the observed SEDs up to far IR wavelengths, but fails to reproduce the radio SED. ... (abstract continues but did not fit here)

Tables

Table 1. Table 1: Details of the IR, sub-mm, and radio observations.

Mission/Telescope	Epoch of observations	${\lambda }_{\text{central}}$/phot. band	Reference
IRAS	1 – 11/1983	12, 25, 60, 100 $\mu$m	Helou & Walker (1988)
VLA	12/1987 – 9/1988, 2/1990, 8/1991	2, 3.6, 6 cm	Taylor et al. (1987, 1990),
			Apparao et al. (1990),
			Dougherty et al. (1991b),
			Dougherty & Taylor (1992)
RAO/NOAO/MKO	4/1988 – 3/1991	$JHKLMN$ bands	Dougherty et al. (1991a)
JCMT	8/1988, 8/1989	0.76, 1.1 mm	Waters et al. (1989, 1991)
IRAM	3/1991, 12/1992	1.2 mm	Wendker et al. (2000)
MSX	1996 – 1997	4.29, 4.35, 8.88, 12.13, 14.65, 21.3 $\mu$m	Egan et al. (2003)
AKARI	5/2006 – 8/2007	9, 18, 65, 90, 140, 160 $\mu$m	Ishihara et al. (2010)
WISE	1/2010 – 11/2010	3.3, 4.6, 12 and 22 $\mu$m	Cutri & et al. (2014)
JVLA	10/2010	0.7, 1.3, 3.5 and 6.0 cm	this work
CARMA	9/2013	3 mm	Klement et al. (2015)
APEX/LABOCA	9/2013, 7–8/2015	0.87 mm	this work

Table 2. Table 2: VLA Calibration Sources.

Source	$\alpha$(J2000.0)	$\delta$(J2000.0)	Date	${\nu }_0$ (MHz)	Flux Density (Jy)
J0102+5824a	01 02 45.7624	58 24 11.137	21-OCT-2010	4.894	1.321$\pm$0.004
			21-OCT-2010	5.022	1.347$\pm$0.004
			30-OCT-2010	8.394	1.965$\pm$0.006
			30-OCT-2010	8.522	1.969$\pm$0.005
			31-OCT-2010	8.394	1.937$\pm$0.003
			31-OCT-2010	8.522	1.948$\pm$0.003
			30-OCT-2010	22.394	2.28$\pm$0.04
			30-OCT-2010	22.522	2.26$\pm$0.03
			31-OCT-2010	22.394	2.08$\pm$0.02
			31-OCT-2010	22.522	2.09$\pm$0.01
			30-OCT-2010	43.342	2.03$\pm$0.10
			30-OCT-2010	43.470	2.04$\pm$0.10
			31-OCT-2010	43.342	1.82$\pm$0.04
			31-OCT-2010	43.470	1.82$\pm$0.04
3C48b	01 37 41.2994	+33 09 35.132	11-OCT-2010	4.895	5.3854∗
			12-OCT-2010	4.895	5.3855∗
			21-OCT-2010	4.895	5.3857∗
			31-OCT-2010	4.895	5.3859∗
			15-OCT-2010	8.395	3.2074∗
			30-OCT-2010	8.395	3.2080∗
			31-OCT-2010	8.395	3.2081∗
			15-OCT-2010	22.395	1.2422∗
			30-OCT-2010	22.395	1.2429∗
			31-OCT-2010	22.395	1.2430∗
			15-OCT-2010	43.215	0.6739∗
			30-OCT-2010	43.343	0.6728∗
			31-OCT-2010	43.343	0.6728∗
J0349+4609c	03 49 18.7416	46 09 59.658	31-OCT-2010	8.394	0.5802$\pm$0.0005
			31-OCT-2010	8.522	0.5772$\pm$0.0005
			31-OCT-2010	22.394	0.410$\pm$0.002
			31-OCT-2010	22.522	0.407$\pm$0.002
			31-OCT-2010	43.342	0.335$\pm$0.005
			31-OCT-2010	43.470	0.333$\pm$0.006
J0359+5057c	03 59 29.7473	50 57 50.162	11-OCT-2010	4.894	8.56$\pm$0.02
			11-OCT-2010	5.022	8.71$\pm$0.02
J0403+2600d	04 03 05.5860	26 00 01.503	31-OCT-2010	4.894	2.344$\pm$0.003
			31-OCT-2010	5.150	2.357$\pm$0.004
			15-OCT-2010	8.394	2.190$\pm$0.003
			15-OCT-2010	8.650	2.187$\pm$0.003
			15-OCT-2010	22.394	1.55$\pm$0.01
			15-OCT-2010	22.650	1.56$\pm$0.01
			15-OCT-2010	43.214	1.15$\pm$0.04
			15-OCT-2010	43.470	1.15$\pm$0.04
J2311+4543e	23 11 47.4090	45 43 56.016	31-OCT-2010	8.394	0.602$\pm$0.001
			31-OCT-2010	8.650	0.605$\pm$0.001
			31-OCT-2010	22.394	0.525$\pm$0.006
			31-OCT-2010	22.650	0.531$\pm$0.007
			31-OCT-2010	43.342	0.43$\pm$0.01
			31-OCT-2010	43.598	0.43$\pm$0.01
J2322+5057e	23 22 25.9822	50 57 51.964	12-OCT-2010	4.894	1.251$\pm$0.008
			12-OCT-2010	4.894	1.245$\pm$0.008

Table 3. Table 3: Deconvolved image characteristics for the new VLA Data.

Source	$\lambda$ [cm]	${\theta }_{\mathrm{FWHM}}$ [arcsec]	PA [degrees]	rms [$\mu$Jy beam-1]
$\gamma$ Cas	6.0	5.4$\times$4.0	$-59.3$	11.0
$\gamma$ Cas∗	3.5	2.9$\times$2.3	$-43.6$	16.0
$\gamma$ Cas∗	1.3	1.1$\times$0.89	$-29.6$	56.4
$\gamma$ Cas∗	0.7	0.63$\times$0.51	$-33.7$	156
$\psi$ Per	6.0	6.7$\times$4.7	+88.7	12.3
$\psi$ Per	3.5	2.6$\times$2.3	+9.8	21.9
$\psi$ Per	1.3	1.1$\times$1.0	+12.4	83.5
$\psi$ Per	0.7	0.67$\times$0.55	+29.7	252
EW Lac	6.0	6.3$\times$4.6	$-89.6$	11.7
EW Lac	3.5	2.7$\times$2.2	+35.1	23.4
EW Lac	1.3	1.2$\times$0.9	+30.6	76.5
EW Lac	0.7	0.60$\times$0.53	+57.0	275
$\eta$ Tau	6.0	4.6$\times$4.1	$-3.2$	9.75
$\eta$ Tau	3.5	2.7$\times$2.2	$-2.5$	23.9
$\eta$ Tau	1.3	1.0$\times$0.9	+14.9	91.1
$\eta$ Tau	0.7	0.62$\times$0.52	+56.9	297.5

Table 4. Table 4: Newly acquired radio observations from JVLA (0.7, 1.3, 3.5 and 6.0 cm).

Star	$F$ at 0.7 cm [mJy]	$F$ at 1.3 cm [mJy]	$F$ at 3.5 cm [mJy]	$F$ at 6.0 cm [mJy]
$\eta$ Tau	2.09 $\pm$ 0.30	0.852 $\pm$ 0.090	0.237 $\pm$ 0.024	0.1432 $\pm$ 0.0098
EW Lac	0.83 (3$\sigma$)	0.351 $\pm$ 0.077	0.133 $\pm$ 0.023	0.089 $\pm$ 0.012
$\gamma$ Cas	1.65 $\pm$ 0.15	0.710 $\pm$ 0.057	0.201 $\pm$ 0.016	0.128 $\pm$ 0.011
$\psi$ Per	3.47 $\pm$ 0.26	1.284 $\pm$ 0.083	0.435 $\pm$ 0.022	0.241 $\pm$ 0.012

Table 5. Table 5: Adopted model parameters.

Star	Spectral typea,b𝑎𝑏a,ba,b𝑎𝑏a,bfootnotemark: $a,b$	Massc𝑐cc𝑐c$c$Secondary gain calibrator for $\psi$ Per. [M⊙]	Distanced𝑑dd𝑑d$d$Secondary gain calibrator for $\eta$ Tau. [pc]	Inclination [°]
$\beta$ CMi	B8 Ve	3.5e𝑒ee𝑒e$e$Secondary gain calibrator for EW Lac.	49.6	43e𝑒ee𝑒e$e$Klement et al. (2015);
$\eta$ Tau	B7 IIIe	3.4	123.6	41f𝑓ff𝑓f$f$Tycner et al. (2005);
EW Lac	B3 IVe-shell	6.1	252	80b,g𝑏𝑔b,gb,g𝑏𝑔b,gfootnotemark: $b,g$
$\psi$ Per	B5 IIIe-shell	4.4	178.9	75hℎhhℎh$h$Delaa et al. (2011);
$\gamma$ Cas	B0.5 IVe	13.2	168.4	42i𝑖ii𝑖i$i$Stee et al. (2012).
$\beta$ Mon A	B4 Ve-shell	5.1	208	70b,g𝑏𝑔b,gb,g𝑏𝑔b,gfootnotemark: $b,g$

Table 6. Table 6: Derived model parameters.

Star	$R_{\text{p}}$ [R⊙]	$R_{\text{e}}$ [R⊙]	$L$ [L⊙]	$\log {\rho }_0$ [g cm-3]	$n$	$R_{\text{out}}$ [R${}_{\text{e}}$]	$R_{\text{c}}$ [R${}_{\text{e}}$]	$E\left(B-V\right)$
$\beta$ CMi	$2.8\pm 0.2$	$4.20$	$185\pm 5$	$-{11.78}_{-0.30}^{+0.18}$	$2.9\pm 0.1$	${40}_{-5}^{+10}$	400	${0.01}_{-0.01}^{+0.02}$
$\eta$ Tau	$8.0\pm 0.5$	$10.56$	$1750\pm 50$	$-{12.23}_{-0.16}^{+0.12}$	$3.0\pm 0.4$	${40}_{-5}^{+10}$	140	$0.06\pm 0.02$
EW Lac	$4.5\pm 0.5$	$5.94$	$3500\pm 100$	$-{11.60}_{-0.18}^{+0.12}$	$2.2\pm 0.1$	${100}_{-20}^{+30}$	360	${0.11}_{-0.03}^{+0.04}$
$\psi$ Per	$5.5\pm 0.5$	$7.26$	$2100\pm 100$	$-{11.48}_{-0.12}^{+0.10}$	$2.5\pm 0.1$	${100}_{-15}^{+5}$	250	${0.08}_{-0.02}^{+0.01}$
$\gamma$ Cas	$7.0\pm 0.5$	$9.24$	$19000\pm 500$	$-{10.48}_{-0.12}^{+0.10}$	$3.3\pm 0.1$	$35\pm 5$	560	${0.04}_{-0.02}^{+0.03}$
$\beta$ Mon A	$4.0\pm 0.5$	$5.28$	$4500\pm 100$	$-{10.60}_{-0.18}^{+0.19}$	${3.0}_{-0.1}^{+0.2}$	$110\pm 40$	380	0.11$\pm$ 0.04

Equations

$$\frac{R_{\text{c}}}{R_{\mathrm{e}}}=\frac{3}{10}{\left(\frac{v_{\text{orb}}}{c_{\text{s}}}\right)}^2,$$

$$\rho (r,z)={\rho }_0{\left(\frac{r}{R_{\text{e}}}\right)}^{-n}\text{exp}\left(-\frac{z^2}{2H^2}\right),$$

$$H(r)=H_0{\left(\frac{r}{R_{\text{e}}}\right)}^{3/2},$$

$${\chi }_{\mathrm{mod}}^2=\sum _{i=1}^N{\left(\frac{\log F_{\text{obs,i}}-\log F_{\text{mod,i}}}{F_{\text{err,i}}/F_{\text{obs,i}}}\right)}^2,$$

Full text

*Aims. *We test the VDD model predictions for the extended parts of a sample of six Be disks that have been observed in the radio to address the question of whether the observed turndown in the spectral energy distribution (SED) can be explained in the framework of the VDD model, including recent theoretical development for truncated Be disks in binary systems.

*Methods. *We combine new multi-wavelength radio observations from the Karl. G. Jansky Very Large Array (JVLA) and Atacama Pathfinder Experiment (APEX) with previously published radio data and archival SED measurements at ultraviolet, visual, and infrared wavelengths. The density structure of the disks, including their outer parts, is constrained by radiative transfer modeling of the observed spectrum using VDD model predictions. In the VDD model we include the presumed effects of possible tidal influence from faint binary companions.

*Results. *For 5 out of 6 studied stars, the observed SED shows strong signs of SED turndown between far-IR and radio wavelengths. A VDD model that extends to large distances closely reproduces the observed SEDs up to far IR wavelengths, but fails to reproduce the radio SED. Using a truncated VDD model improves the fit, leading to a successful explanation of the SED turndown observed for the stars in our sample. The slope of the observed SEDs in the radio is however not well reproduced by disks that are simply cut off at a certain distance. Rather, some matter seems to extend beyond the truncation radius, where it still contributes to the observed SEDs, making the spectral slope in the radio shallower. This finding is in agreement with our current understanding of binary truncation from hydrodynamical simulations, in which the disk does extend past the truncation radius. Therefore, the most probable cause for the SED turndown is the presence of binary companions that remain undetected for most of our sources.

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 ($\lesssim 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 $\sim$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 $\mu$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 $\rho\propto 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, $\zeta$ Tauri (Carciofi et al., 2009), $\beta$ 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_{\text{c}}$ was derived by Krtička et al. (2011) for isothermal disks:

[TABLE]

where $R_{\mathrm{e}}$ is the stellar equatorial radius, $v_{\text{orb}}$ is the disk orbital speed above the stellar equator and $c_{\text{s}}$ is the speed of sound. The typical values of $R_{\text{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., submitted).

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: $\varphi$ 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 ($\sim 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., $\gamma$ Cas), the H$\alpha$ 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.

\onltab

\onltab

3.1 Ultraviolet spectra

For the UV part of the spectrum, data from the International Ultraviolet Explorer (IUE)333https://archive.stsci.edu/iue/ and the Wisconsin Ultraviolet Photo-Polarimeter Experiment (WUPPE)444https://archive.stsci.edu/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 Analyser555http://svo2.cab.inta-csic.es/theory/vosa/ (VOSA; Bayo et al., 2008) and the VO Spectral Analysis Tool666http://esavo.esa.int/vospec/ (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 $\beta$ 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 $\beta$ 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 $\beta$ Mon A, we used a new APEX measurement from the bolometer camera LABOCA (Siringo et al., 2009) at the wavelength of 870 $\mu$m. The observed flux is $23.5\pm 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\pm 3$ arcsec. Although close to the edge of the FWHM, the B component may still slightly contribute to the observed flux.

For $\beta$ CMi, we also have an APEX/LABOCA measurement, that was originally published in Klement et al. (2015). However, we newly reduced the data using Crush777http://www.submm.caltech.edu/~sharc/crush/index.html version 2.32-1 (Kovács, 2008), with the updated flux value for $\beta$ CMi being $38.6\pm 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 $\beta$ 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 4). 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 $\sim 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 $\sim$65, 13, 4, and 3.5 minutes at 6.0 cm, 3.5 cm, 1.3 cm, and 0.7 cm, respectively, with the exception of $\gamma$ 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 3C48 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 Cookbook888http://www.aips.nrao.edu/cook.html. 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 ($\sim$1 minute) 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 3C48 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 $\sim$10 seconds. 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 2 (online only).

Imaging was performed using the AIPS task IMAGR with robust +1 weighting of the visibilities. Table 3 (online only) 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 CADC999http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/. ESPaDOnS is mounted at the CFHT on Hawaii and described by Donati et al. (2006).

• •

The TBL/Narval archive101010tblegacy.bagn.obs-mip.fr/narval.html. 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 spectropolarimeter111111http://archive.stsci.edu/hpol/ 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_{\text{eff}}$ and $\log{g}$.

The rotating central star is thus completely described by four independent parameters: mass $M$, polar radius $R_{\text{p}}$, luminosity $L$, and rotation rate $W$. The $W$ parameter is defined as $W=v_{\text{rot}}/v_{\text{orb}}$, where $v_{\text{rot}}$ is the rotational velocity at the stellar surface and $v_{\text{orb}}$ is the Keplerian orbital velocity directly above the stellar surface. The gravity darkening exponent $\beta$ that appears in the classical form of the gravity darkening law, $T_{\text{eff}}\propto g^{\beta}$ (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:

[TABLE]

where $r$ and $z$ are respectively the radial and vertical cylindrical coordinates, $\rho_{0}$ is the density of the disk close to the stellar surface, $R_{\text{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:

[TABLE]

where $H_{0}=c_{\text{s}}v_{\text{orb}}^{-1}R_{\text{e}}$ is the pressure-supported scale height at the base of the disk, $c_{\text{s}}=[(k_{\text{B}}T_{\text{k}})/(\mu m_{\text{a}})]^{1/2}$ is the isothermal sound speed, $k_{\text{B}}$ is the Boltzmann constant, $T_{\text{k}}$ is the gas kinetic temperature, $\mu$ is the mean molecular weight, and $m_{\text{a}}$ is the atomic mass unit. $T_{\text{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\approx 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_{\text{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_{\text{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. (2011).

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 ($\rho_{0}$, $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_{\text{e}}/R_{\text{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 5.

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 module141414available online at http://dan.iel.fm/emcee under the MIT License (Foreman-Mackey et al., 2013).

We proceed to the remaining physical parameters of the central star, namely the polar radius $R_{\text{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_{\text{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_{\text{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_{\text{p}}$ and $L$ parameter combination was determined via $\chi^{2}$ minimization, where for the best-fit model $\chi^{2}=\chi^{2}_{\text{min}}$. Since we were dealing with SED modeling, where both the wavelengths and fluxes span many orders of magnitude, we used the logarithmic version of $\chi^{2}$:

[TABLE]

where $F_{\text{obs}}$ is the observed flux, $F_{\text{mod}}$ is the model flux, $F_{\text{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 $\Delta\chi^{2}_{0.683}=\chi^{2}_{\mathrm{mod}}-\chi^{2}_{\mathrm{min}}$ contour corresponding to the $68.3\%$ confidence level (1-$\sigma$), where $\Delta\chi^{2}_{0.683}$ is a function of the number of degrees of freedom of the fit (we refer to, e.g., Chapter 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 $\rho_{0}$, the density exponent $n$, and the disk outer radius $R_{\text{out}}$. We set an upper limit to the $R_{\text{out}}$ parameter, so that it cannot exceed the value of the critical radius $R_{\text{c}}$ estimated by means of Eq. (1), in which the sound speed $c_{\text{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_{\text{out}}$ to an intermediate value of 50 stellar radii, and a grid of models with $\rho_{0}$ and $n$ in the range 1$\times 10^{-13}$ – 1$\times 10^{-10}$ g cm*-3* and 1.5 – 5.0, respectively, was computed. To determine the parameters $\rho_{0}$ 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 $\rho_{0}$ 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 $\rho_{0}$ and $n$ were determined in the same way as for the central star parameters $R_{\text{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_{\text{p}}$ and $L$ in a straightforward way. We adopted an iterative modeling procedure for such stars: A first look estimate for $R_{\text{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 $\rho_{0}$ 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 $\rho_{0}$ and $n$, we can proceed to the determination of the $R_{\text{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 $\gtrsim 20$ $R_{\text{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., $\beta$ CMi). In the latter case, the sub-mm/mm region was included when searching for the best-fit values of $\rho_{0}$ and $n$. Measurements corresponding to upper limits were used for visual checking, but not for the modeling procedure itself. The uncertainty of the $R_{\text{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 $\beta$ CMi from Klement et al. (2015). Then we present the results for stars for which we have new JVLA measurements ($\eta$ Tau, EW Lac, $\psi$ Per, $\gamma$ Cas), followed by the results for $\beta$ 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 6. In Table 6 we also include the values of $R_{\text{c}}$ calculated using Eq. (1).

5.1 $\beta$ CMi (HD 58715; HR 2845)

The star $\beta$ 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 $\rho_{0}$, $n$ and $R_{\text{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 $\rho_{0}$ 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 $35^{+10}_{-5}$ $R_{\text{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 $40^{+10}_{-5}$ $R_{\text{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 $\eta$ 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\pm 2.1$ km/s. Heros measurements around MJD 52000 are on average approximately 4 km/s lower than ESPaDOnS measurements at MJD 55000, but given the different instrument parameters, the values are well compatible and within 3-$\sigma$ 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$\alpha$ emission was at a constant level of $I/I_{\mathrm{c}}=2.2$, and $W_{\lambda}=-3.5$ Å until approximately 2005, and since then steadily decreased by a small but detectable amount to $I/I_{\mathrm{c}}=2.0$, and $W_{\lambda}=-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-$\sigma$ of one another.

The results of the SED modeling of $\eta$ Tau are plotted in Fig. 2. The non-shell nature of $\eta$ Tau allows us to use the UV spectrum to determine the central star parameters $R_{\text{p}}$ and $L$ without any complications. Moreover, the disk of $\eta$ Tau is apparently very tenuous and starts to contribute to the observed SED only at around 5 $\mu$m. We thus use both the UV spectrum and the visual photometry to constrain $R_{\text{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 $\mu$m to mm was used to find the best-fit values of the disk parameters $\rho_{0}$ 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 $R_{\text{out}}=40^{+10}_{-5}$ $R_{\text{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 $\geq$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 51000 to 52500, are at approximately $-22$ km/s, while ESPaDOnS values, obtained around MJD 55400, are at $-28$ km/s. 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$\alpha$ 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 $\eta$ 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 $\sigma$ 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 $\mu$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 $\rho_{O}$ and $n$, while the mm and radio part, sensitive to varying the disk size, was used to determine $R_{\text{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 $100^{+30}_{-20}$ $R_{\text{e}}$. Moreover, the new JVLA data set reveals a similar result as given above for $\eta$ 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 $\psi$ Per (HD 22192; HR 1087)

The star $\psi$ 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$° $\pm$ $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. While the RV is varying by approximately 1500 m/s 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 $\psi$ 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_{\mathrm{c}}=6$ and $W_{\lambda}=-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 $\psi$ Per are plotted in Fig. 4. Although $\psi$ 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_{\text{p}}$ and $L$. The disk contribution to the observed SED starts to be significant at near-IR wavelengths. Although the value of $\rho_{0}$ 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 $100^{+5}_{-15}$ $R_{\text{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, $\psi$ Per is the only star that has been angularly resolved in radio ($\lambda$ = 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_{\text{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${}_{\text{out}}=100$ $R_{\text{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 $\mu$Jy) is relatively close to the measured value ($804\pm 60$ $\mu$Jy). The model with the disk size of $500$ $R_{\text{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\times 20.1$ mas and $85.2\times 47.8$ mas for the best-fit and the $500$ $R_{\text{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.

5.5 $\gamma$ Cas (HD 5394; HR 264)

The star $\gamma$ 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). $\gamma$ 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\pm 0.09$ km/s and a period of $203.52\pm 0.08$ days. $\gamma$ 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).

$\gamma$ 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), $\gamma$ 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$\alpha$ line also seems to indicate binarity, since it often shows no clear peaks or symmetry, which are signs of disturbances in the disk.

Judging from the H$\beta$ profiles, which behave closer to the expectations for a classical Be star than the H$\alpha$ profiles, $\gamma$ 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 $\gamma$ 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_{\mathrm{c}}$ of H$\alpha$ varied, but remained in the range between approximately 4 and 5. In other words, while the disk of $\gamma$ Cas is not stable in the same sense as the one of, for instance, $\psi$ 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 $\sigma$ below those of AKARI/WISE. The same is true for the observed mm fluxes, as the JCMT and IRAM measurements are discrepant by several $\sigma$. This suggests mild variability in the disk, consistent with what is observed in the H$\alpha$ 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 $\gamma$ 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_{\text{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 $\rho_{0}$ 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 $35^{+5}_{-5}$ $R_{\text{e}}$. However, the observed slope of the radio SED is again flatter than the slope of the truncated model.

5.6 $\beta$ Mon A (HD 45725; HR 2356)

$\beta$ 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., 2016); 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 51000 and 56000. The narrow circumstellar line Fe II 5169 shows small RV variations with an amplitude of approximately 3 km/s. Again, however, it is not clear whether this variation is due to binarity or changes in the disk, since $\beta$ Mon A is a V/R variable Be star.

No HPOL data are available for $\beta$ Mon A. Judging from the line emission, the disk is generally strong, with values of $I/I_{\mathrm{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 $\gamma$ 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 $\sigma$ below the detection, suggesting a mild variability of the radio fluxes.

The results of the SED modeling of $\beta$ Mon A are plotted in Fig. 7. $\beta$ 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 $\beta$ 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 $\beta$ 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-$\sigma$ 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\pm 40\,R_{\mathrm{e}}$, which is significantly lower that the critical radius value of $260\,R_{\mathrm{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 $\rho_{0}$ 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\geq 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\lesssim 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{\rho_{0}}$ 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 $\psi$ 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 $\beta$ 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.

We found that the disk truncation occurs at much smaller distances from the central star than where the supersonic point (the critical radius, $R_{\text{c}}$, in Table 6) 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_{\mathrm{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_{\mathrm{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 $\beta$ 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$\alpha$ line led to the detection of the binary, with an orbital period of 170 days and RV amplitude of 2.25 km/s (Dulaney et al., 2017). The value of the detected period, if a circular, coplanar orbit and a 1 $M_{\odot}$ companion are assumed, indicates a semi-major axis of $\sim 50$ $R_{\text{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 $\sim 25$ $R_{\text{e}}$. The result that the disk is truncated at $\sim 40$ $R_{\text{e}}$ rather favors the 3:2 resonance with the orbit as the location of the truncation radius for the case of $\beta$ CMi.

For the other confirmed binary in our sample, $\gamma$ Cas, the semi-major axis of the orbit is $a\simeq 1.64$ au (Nemravová et al., 2012), corresponding to $38\,R_{\mathrm{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 ($\eta$ Tau, EW Lac, and $\psi$ 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 derived for $\beta$ 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 the observed SEDs covering the wavelength interval from the UV to the far-IR. This further establishes the VDD model as the best model 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\sim 3$, typically associated with disks in steady-state) than shallower profiles ($n\sim 2$, usually found in Be stars whose disks are dissipating). However, for these latter cases (namely EW Lac and $\psi$ 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 $\beta$ Mon A, for which the measurements are not conclusive. The effect is present both for confirmed binaries ($\beta$ CMi, $\gamma$ 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, $\psi$ 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_{\text{e}}$ that shows good agreement with the observed shape of the visibility curve. However, this model clearly overestimates the observed flux density.

Acknowledgements.

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).