A Multi–Observing Technique Study of the Dynamical Evolution of the Viscous Disk around the Be Star ω CMa

ω (28) CMa (HD 56139, HR 2749; HIP 35037; B2 IV-Ve) is one of the brightest Be stars in the sky (with m_v ≈ 3.6 to 4.2 mag), and it has caught the attention of observers for more than five decades. It is a nearly pole-on star (i ≈ 15°) so the measured projected rotational velocity of 80 km s⁻¹ (Slettebak et al. 1975) is only a fraction of the true equatorial velocity, estimated to be 350 km s⁻¹ (Maintz et al. 2003). A 1.37 day line-profile variability has been observed in ω CMa, suggesting that it is a nonradial pulsator (Baade 1982). Later, this was confirmed by Štefl et al. (1999) and Maintz et al. (2003) by studying the line-profile variations caused by NRP for various photospheric absorption lines of different species, including Balmer lines, He i, Mg ii, and Fe ii. Recently, from space photometry with BRITE-Constellation (Weiss et al. 2014), Baade et al. (2017) found that the 0.73 day⁻¹ frequency (corresponding to the 1.37 day period) is part of an NRP frequency group between ∼0.55 day⁻¹ and ∼0.8 day⁻¹. Another frequency group between ∼1.15 day⁻¹ and ∼1.45 day⁻¹ seemed to exhibit a much increased amplitude at a time when the mean brightness also increased, i.e., matter was ejected into the disk when the NRP amplitude was high. Observations with the Transiting Exoplanet Survey Satellite of hundreds of Be stars have established such a correlation in dozens of other Be stars (Labadie-Bartz et al. 2020), suggesting that nonlinear coupling of NRP modes (Baade et al. 2018b) can indeed lead to mass ejection events in Be stars. The stellar parameters of ω CMa used in this study are summarized in Table 1 and were obtained by Maintz et al. (2003).

Table 1. Stellar, Disk, and Geometrical Parameters of ω CMa

	Parameter	Value	Reference
Input Parameters for Modeling
Star	M	9.0 M⊙	Maintz et al. (2003)
	L	5224 L⊙	Maintz et al. (2003)
	Tpole	22,000 K	Maintz et al. (2003)
	Rpole	6.0 R⊙	Maintz et al. (2003)
	Req	7.5 R⊙	Maintz et al. (2003)
	log gpole	3.84	Maintz et al. (2003)
	vrot	350 km s−1	Maintz et al. (2003)
	W	0.73	Maintz et al. (2003)
Disk	${\dot{M}}_{\mathrm{inj}}$ (min)	2.0 × 10−10 M⊙ yr−1	G18
	${\dot{M}}_{\mathrm{inj}}$ (max)	3.7 × 10−7 M⊙ yr−1	G18
	$-{\dot{J}}_{* ,\mathrm{std}}$ (min)	2.9 × 1033 g cm2 s−2	G18
	$-{\dot{J}}_{* ,\mathrm{std}}$ (max)	5.4 × 1036 g cm2 s−2	G18
	T0	13,200 K	Carciofi et al. (2012)
	Rout	1000 Req	G18
Other parameters
	i	12°–18°	this work
	vcrit	436 km s−1	Maintz et al. (2003)
	${m}_{{\rm{v}}}^{* }$	4.22 ± 0.05	G18
	d	${280}_{-11}^{+13}$ pc	Gaia Collaboration et al. (2016, 2020)

Note. ${m}_{{\rm{v}}}^{* }$ is the magnitude of the star during the diskless phase.

Download table as:  ASCIITypeset image

Carciofi et al. (2012) used SINGLEBE and HDUST to study the disk dissipation of ω CMa between 2003 and 2008. With time-dependent models of the dissipating disk they determined the α parameter. Moreover, they showed that the stellar wind is not a probable mechanism for AM injection into the disk. Later, their work was followed up by Ghoreyshi et al. (2018, hereafter G18), who presented a model of the full V-band light curve, spanning more than 30 yr of data. The model addressed any photometric variability larger than 0.05 mag and longer than 2 months. Any shorter variability was excluded because the contribution to the total gas content of the disk is small and is usually poorly sampled.

G18 showed that the VDD model could reproduce the data well. It was determined that α changes during the different epochs of the disk life as $\dot{J}$ varies over time. They also found that the light curve could only be reproduced if quiescence phases are interpreted as a reduction of $\dot{J}$, rather than a complete cessation of it, as is commonly assumed (see Carciofi et al. 2012).

In G18, V-band photometric data were investigated in detail, and the results relevant to this study are summarized here. Since 1982, ω CMa has exhibited quasi-regular cycles, each lasting between about 7.0 and 10.5 yr. Each cycle consists of two main parts: (1) an outburst phase represented by a fast increase in brightness and a subsequent plateau (this increase is not always smooth, and the peak brightness plateau lasts about 2.5–4.0 yr) and (2) a quiescence phase lasting about 4.5–6.5 yr that is characterized by a fast decline in brightness and a subsequent slow fading. During these phases the brightness of the system in the V band changes from about 03 to 05. Throughout this text we refer to the cycles as Ci and to the phases as Oi and Qi for outburst and quiescence, ⁶ respectively, where i is the cycle number. All four cycles are labeled in Figure 1.

Zoom In Zoom Out Reset image size

Because the V-band excess comes from the innermost part of the disk (Carciofi 2011), the numerical solutions used in G18 were not constrained past a few stellar radii from the star. In this paper, the same parameters of the model of G18 are used to study other observables of ω CMa, including polarimetric, spectroscopic, and photometric data at other wavelengths. Our goal is to compare the predictions of the G18 model for a much larger volume of the disk out to 30 stellar radii. For this work, and following G18, only variability longer than 2 months is investigated.

Different observables emanate from different parts of the disk (Carciofi 2011). For example, continuum polarization originates near the star, and spectral lines form in various locations based on the wavelength-dependent opacity and the source function for the line. Therefore, probing Be star disks with a variety of observational techniques and in a variety of wavelength regimes allows us to perform comprehensive studies of these systems.

The main goal of this paper is to use the VDD model to study the temporal variations seen in ω CMa's data by a variety of observational techniques. Reproducing the observed data in a wide range of wavelengths is a new challenge for the VDD model that has not been performed.

Section 2 describes the observational data available for ω CMa. Section 3 presents the theoretical concepts that are used in this work. In Section 4 the modeling of the available multi–observing technique data is presented. In Section 5 possible solutions to the discrepancies between the data and models are discussed. In Section 6 our conclusions and plans for future work are presented.

At the end of 2008, Sebastian Otero alerted the community (in a private communication to our now-deceased colleague Stanislav Štefl) to a new outburst that had begun; thus a broad suite of observations was undertaken. In addition to visual photometry, JHKL photometry was obtained with the Mk II photometer (Glass 1973) of the South African Astronomical Observatory (SAAO) and the CAIN-II (Camara Infrarroja) Tenerife/Telescopio Carlos Sánchez camera (Cabrera-Lavers et al. 2006), and Q1- and Q3-band measurements were made with VISIR (Lagage et al. 2004) on the Very Large Telescope (VLT) at the European Southern Observatory (ESO). Finally, we also included photometric data in the UBV Johnson and uvby Strömgren (Strömgren 1956; Crawford 1958) filters collected in the Long-term Photometry of Variables (LTPV) project (Manfroid et al. 1995) during C1.

The above campaign also produced optical echelle spectra from the Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) and VLT (2008 October–2009 March), the Fiber-fed Extended Range Optical Spectrograph (feros) at La Silla (Kaufer et al. 1999), and the 1.6 m telescope at Observatório Pico dos Dias (OPD; 2009 January–20164ri) initially using the ECASS spectrograph ⁷ and since 2012 using the MUSICOS spectrograph (Baudrand & Böhm 1992).

In addition to the observational effort described above for C4 (see Figure 1), we obtained data for some of the previous cycles in the literature. For C1 and C2, we acquired spectroscopy from the Short-wavelength Prime camera of the International Ultraviolet Explorer (IUE) ⁸ and the Heidelberg Extended Range Optical Spectrograph (heros) ⁹ /feros, respectively. For C3, we found spectroscopy from the Coudé Echelle Spectrometer (CES), ¹⁰ feros, the Lhires spectroscope in Observatoire Paysages du Pilat, ¹¹ and Ondrejov Observatory. ¹² Additional spectroscopic data for C4 came from the Echelle Spectropolarimetric Device for the Observation of Stars (ESPaDOnS; Donati 2003), OPD, PHOENIX (Hinkle et al. 1998), Ritter Observatory, ¹³ and UVES.

Figure 3 provides an example of the observed hydrogen lines of ω CMa. The top panel shows the V-band photometric data with the date the spectroscopic data were observed, indicated by the colored vertical solid line. The bottom panels show the flux relative to the local continuum for the four main hydrogen lines. Usually the peak-emission-to-continuum ratio (E/C) of the Hα and Hβ lines is highest at the end of quiescence and lowest during the outburst. This seemingly contradictory behavior is well explained by the models as will be seen in Section 4.3.

Zoom In Zoom Out Reset image size

BVRI imaging polarimetry was obtained with the 0.6 m telescope at OPD (Magalhães et al. 2006). Reduction of the OPD polarimetric data, observed by the IAGPOL instrument, followed the standard procedures outlined by Magalhães et al. (1984, 1996) and Carciofi et al. (2007). IAGPOL has an instrumental polarization lower than about 0.005% (Carciofi et al. 2007). The middle panel of Figure 4 shows the polarimetric data of C4 of ω CMa in the B, V, R, and I filters, alongside the V-band photometric data (top panel) and polarization angle, θ, measured east from celestial north (bottom panel). The polarization level is very small, as expected from the fact that ω CMa is observed at rather small inclination angles (Halonen & Jones 2013). This happens because at nearly pole-on orientations an axisymmetric disk will appear as an almost circular structure, which results in the cancellation of the polarization vectors. The small observed polarization level indicates that the interstellar (IS) polarization is likely to be very low as well. Both combined factors make analysis of the polarization data uncertain, as will be discussed in Section 4.1. The observational logs of the spectroscopic and polarimetric data used in this paper are listed in Appendix A and in Tables 3 and 4, respectively. Also, the data are available at CDS. ¹⁴

Zoom In Zoom Out Reset image size

SINGLEBE solves the isothermal 1D time-dependent fluid equations (Pringle 1981) in the thin-disk approximation and provides the disk surface density, Σ(r, t).

The 1D grid used in the SINGLEBE code models the disk between R_eq (the equatorial radius of the star) and R_out (the outer radius of the disk). The grid is a logarithmic array with an optional number of cells. One cell is arbitrarily selected as the location where mass from the central star is entirely injected, ${\bar{r}}_{\mathrm{inj}}{R}_{\mathrm{eq}}$, where ${\bar{r}}_{\mathrm{inj}}$ is a dimensionless quantity. The viscosity parameter, as a function of time, α(t), and ${\dot{J}}_{* ,\mathrm{std}}(t)$ are the input parameters of the code. SINGLEBE determines how the injected matter spreads in the disk. More details about SINGLEBE can be obtained in the original publication (Okazaki 2007), and a description of the boundary conditions adopted can be found in Rímulo et al. (2018).

The Monte Carlo radiative transfer code HDUST is a fully three-dimensional code that simultaneously solves the radiative equilibrium, radiative transfer, and non-LTE statistical equilibrium equations to obtain the ionization fraction, hydrogen-level populations, and electron temperature as a function of position in a 3D envelope around the star (Carciofi et al. 2004; Carciofi & Bjorkman 2006, 2008). In order to convert the surface density provided by SINGLEBE to volume density, HDUST uses a Gaussian vertical density profile with a 1.5 power-law isothermal disk scale height. With these quantities, HDUST produces the emergent spectral energy distribution, including the emission line profiles, as well as the polarized spectrum and synthetic images.

To date, HDUST has been used in a variety of theoretical studies of Be stars' disks (e.g., Carciofi & Bjorkman 2006, 2008; Haubois et al. 2012, 2014; Faes et al. 2013). The most relevant study for this work is that by Haubois et al. (2012), who used hydrodynamic simulations to show the capability of the VDD model for reproducing the light curves of Be stars. In addition, HDUST has been used in several studies where the model predictions are constrained by observations such as visible and IR photometry (e.g., Carciofi et al. 2012; Baade et al. 2018a; Ghoreyshi et al. 2018; Rímulo et al. 2018), radio photometry (e.g., Klement et al. 2017, 2019), polarimetry (e.g., Carciofi et al. 2007, 2009; Faes et al. 2016), spectroscopy (e.g., Carciofi et al. 2010; Suffak et al. 2020), and spectro-interferometry (e.g., Carciofi et al. 2009; Klement et al. 2015; Faes et al. 2016; de Almeida et al. 2020).

As mentioned earlier, the central star is oblate. Consequently, the polar regions of the star have a greater effective gravity than the equatorial regions. According to the von Zeipel (1924) theorem, this latitudinal dependence of the effective gravity causes a latitudinal dependence of the flux, in the sense that the poles are brighter (hotter) and the equator darker (cooler). This gravity darkening effect plays a key role in determining the surface distribution of the flux and therefore the latitudinal dependence of the temperature.

In its original formalism, applicable for a purely radiative envelope, the von Zeipel theorem can be written as

$$\begin{eqnarray}&&{T}_{\mathrm{eff}}(\theta )\propto {g}_{\mathrm{eff}}^{\beta }(\theta ),\end{eqnarray} \tag{ 3 }$$

with β = 0.25. However, interferometric studies of stellar spectral classes F to B have suggested that the β parameter is typically in the range of 0.18–0.25 (van Belle et al. 2006; Monnier et al. 2007; Che et al. 2011; van Belle 2012; Domiciano de Souza et al. 2018; Hadjara et al. 2018) with the most likely value being 0.21 (van Belle 2012). The theoretical study of Espinosa Lara & Rieutord (2011) suggested that the value of β is a function of the rotational rate of the star. Following these authors, G18 adopted a β of 0.19 for ω CMa. The same value is used in this paper. Using the β parameter, the polar radius of the star (R_pole), and the critical fraction of the rotational velocity (W; see Rivinius et al. 2003 for how this parameter is defined) as input parameters, HDUST calculates the geometrical oblateness and gravitational darkening of the star.

Although HDUST has the ability to take into account the opacity of dust grains (e.g., for B[e] stars, see Carciofi et al. 2010), our models were calculated for dust-free gaseous disks, consistent with our current understanding of Be stars.

The linear polarization level can be expressed in terms of the Stokes parameters, Q and U (Clarke 2010), as

$$\begin{eqnarray}&&P=\sqrt{{Q}^{2}+{U}^{2}}.\end{eqnarray} \tag{ 4 }$$

The polarization position angle is

$$\begin{eqnarray}&&\theta =\displaystyle \frac{1}{2}\arctan \left(\displaystyle \frac{U}{Q}\right).\end{eqnarray} \tag{ 5 }$$

One common issue regarding the interpretation of polarimetric data is the removal of the IS contribution to the observed signal. The observed polarization, decomposed in its Stokes Q and U parameters, can be written as

$$\begin{eqnarray}&&{Q}_{\mathrm{obs}}={Q}_{\mathrm{IS}}+{Q}_{\mathrm{int}},\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{U}_{\mathrm{obs}}={U}_{\mathrm{IS}}+{U}_{\mathrm{int}},\end{eqnarray} \tag{ 7 }$$

i.e., without knowing the IS components (Q_IS and U_IS) of the observed polarization, the components (Q_int and U_int) of intrinsic polarization are unknown. This is shown schematically in the top panel of Figure 5. Measuring the IS component of the polarization can be a challenging task (e.g., Wisniewski et al. 2010). As there is no reliable information about this quantity for ω CMa in the literature, we employ below three different methods to determine the position angle of the intrinsic polarization and to estimate the IS polarization (θ_IS and P_IS).

Zoom In Zoom Out Reset image size

4.1.1.  Q − U Method

The first method explores the fact that the intrinsic polarization is variable, while on the same timescale and at a given wavelength, P_IS is not. We begin examining the bottom panel of Figure 5, which shows, in a schematic way, how the process of formation and dissipation of a Be disk appears in the Q − U diagram. The intrinsic polarization angle on the sky is θ_int. When P_int is zero (no disk), the observed polarization will be due solely to the IS component (P_obs = P_IS). As the disk grows (and dissipates), the magnitude of P_int changes, but not the angle (assuming that the disk is axisymmetric and lies along the equatorial plane). This is shown in the bottom panel of Figure 5 as a track of points along the 2θ_int direction, which indicates that the angle of the track is a measure of θ_int (Draper et al. 2014). More specifically, θ_int should be parallel to the minor elongation axis of the Be disk (we note that the disk may not be elliptic but appears as an ellipse in the plane of the sky in the line of sight of the observer if it is not seen pole-on). In the case of disks confined to the equatorial plane, θ_int also describes the position angle of the spin axis of the star, measured east from celestial north.

Figure 6 shows the Q − U diagram of the polarization data of C4. The original data have some individual points with significant variations and large error. Therefore, we binned the data in time intervals of 100 days. For all four bands, the measurements form a straight path in the Q − U diagram as explained above. A simple linear least-squares regression fit (solid red lines in Figure 6) indicates that the angle of this path is 117° ± 76, 102° ± 58, 110° ± 46, and 104° ± 48 for B, V, R, and I, respectively, which means that θ_int should be half of these values. The errors were estimated using ±1σ uncertainty. These numbers are listed in Table 2 with a QU superscript, to indicate that they were obtained using the Q − U diagram method. The average value of the intrinsic position angles for the four filters is 542 ± 29, where the uncertainty is the standard deviation of the mean.

Zoom In Zoom Out Reset image size

Table 2. Different Estimates of the IS Polarization and the Position Angle of the Intrinsic Polarization

Method	Parameter	B Filter	V Filter	R Filter	I Filter	Average
QU	${\theta }_{\mathrm{int}}^{{QU}}$ (°)	58.4 ± 3.8	51.2 ± 2.9	55.2 ± 2.3	52.0 ± 2.4	54.2 ± 2.9
	${P}_{\mathrm{IS}}^{\mathrm{LC}}$ (%)	0.12 ± 0.01	0.15 ± 0.01	0.11 ± 0.01	0.11 ± 0.01	0.12 ± 0.02
LC	${\theta }_{\mathrm{IS}}^{\mathrm{LC}}$ (°)	56.4 ± 1.5	53.5 ± 1.8	56.0 ± 2.1	56.6 ± 1.9	55.6 ± 1.2
	${\bar{\theta }}_{\mathrm{int}}^{\mathrm{LC}}$ (°)	56.2 ± 10.9	62.9 ± 9.8	63.1 ± 2.8	58.8 ± 9.3	60.2 ± 2.9
	${P}_{\mathrm{IS}}^{\mathrm{field}}$ (%)	0.13 ± 0.02	0.12 ± 0.03	0.11 ± 0.01	0.08 ± 0.01	0.11 ± 0.02
Field	${\theta }_{\mathrm{IS}}^{\mathrm{field}}$ (°)	59.7 ± 4.6	53.5 ± 6.6	57.9 ± 2.9	60.8 ± 1.6	58.0 ± 2.8
	${\bar{\theta }}_{\mathrm{int}}^{\mathrm{field}}$ (°)	54.5 ± 11.3	62.7 ± 10.7	58.2 ± 8.0	55.0 ± 11.6	57.6 ± 3.3

Notes. QU, LC, and Field indicate the Q − U diagram, the light curve, and the field star, respectively, as the methods that have been used to estimate the intrinsic polarization of ω CMa. The numbers in bold give the average value of the intrinsic position angles of four different filters measured by each particular method.

Download table as:  ASCIITypeset image

The Q − U method requires that the position angle of the disk remain relatively constant over time. We can estimate the validity of this assumption by measuring the correlation between the Stokes Q and U parameters. ρ ≈ −1 as measured, e.g., by the Pearson correlation coefficient method indicates a strong correlation. A weak correlation, i.e., ρ ≈ 0, means the disk behavior is complex and the validity of the method is compromised. We obtained ρ = −0.34, −0.51, −0.35, and −0.60 for the B, V, R, and I filters, respectively, which imply that the errors derived by the linear least-squares regression fit are underestimated. However, our results suggest that the intrinsic polarization angle of the star might be close to 52°, corresponding to the correlation coefficient closest to −1 (I band).

4.1.2. Light-curve Method

The second method uses the light curve itself. The results of G18 indicate that at the end of C4, the V-band excess is very small. According to the best G18 model fits, the inner disk at that phase is very tenuous (Figure 1; see also second panel of Figure 11 of G18). Therefore, if one assumes that P_int at that phase is very low, the observed polarization should be very close to P_IS. According to the relative proximity of ω CMa, small values for P_IS are expected (e.g., Serkowski et al. 1975; Yudin 2001).

Figure 7 presents the photometric and polarimetric data of ω CMa in the V band. The observed polarization is shown as gray circles. The two vertical red lines in Figure 7 indicate the boundaries of the phase assumed for the star to be almost diskless. An average of the data at this phase (a single point for the B, V, and R filters and a few points for the I filter; see Figure 18 in Appendix B) gives us ${Q}_{\mathrm{IS}}^{\mathrm{LC}}$ and ${U}_{\mathrm{IS}}^{\mathrm{LC}}$, which in turn provide the values for ${P}_{\mathrm{IS}}^{\mathrm{LC}}$ and ${\theta }_{\mathrm{IS}}^{\mathrm{LC}}$, listed in Table 2 for each filter. Here the superscript LC indicates that the estimates were made using the light curve itself. Using Equations (4)–(7) and the estimated value of IS polarization, the intrinsic polarization of ω CMa, P_int and θ_int, was calculated. The average value of ${\bar{\theta }}_{\mathrm{int}}^{\mathrm{LC}}$ is 602 ± 29, in good agreement with the value estimated using the Q − U method. We defer for later a discussion on the intrinsic polarization levels.

Zoom In Zoom Out Reset image size

4.1.3. Field Star Method

Finally, we estimated P_IS using the field star method, by which one or more stars that are physically close to the target star and that are known to have no intrinsic polarization are used as a proxy of the IS polarization. By IAGPOL, we observed HD 56876, a B5Vn star (Houk 1982) with m_v ≈ 6.4, a Gaia distance of ${282}_{-3}^{+4}$ pc, and an angular distance of 067 from ω CMa. The distance of ω CMa inferred from the parallax measured by Gaia is ${280}_{-11}^{+13}$ pc (Gaia Collaboration et al. 2016, 2020).

The results of this method are listed in Table 2 with a "field" superscript. Agreement with the previous method based on the light curve is quite good, as both results have very similar values for P_IS and θ_IS. The three estimates of θ_int (shown as bold numbers in Table 2) also agree. By modeling the Brγ interferometric data for ω CMa, Štefl et al. (2011) showed that the position angle of the major axis of the star's disk is −29° (see the lower left panel of Figure 1 in their paper). Since there is a 90° difference between the elongation of the minor and major axes, the position angle for the minor axis of the disk is 61°, which is in agreement with the results presented here.

The estimated P_int level from the field star method is illustrated with blue circles in the bottom panel of Figure 7. A positive correlation between the polarization level and the brightness of the star can be seen. It appears that the polarization level follows the variation of the V-band photometric data with a lag; for instance, in the dissipation phase the drop in polarization is slower than the drop in brightness. This agrees with the theoretical studies of Haubois et al. (2014). This behavior is easily explained when we consider that the V-band polarization probes a slightly larger volume of the disk than the V band (see Figure 1 of Carciofi 2011) and therefore, the timescales for viscous dissipation will be longer. There is also some intrinsic scatter in the data, which is likely the result of disk variability. As shown by Haubois et al. (2014), the polarization level responds quickly to changes in $\dot{J}$. The flickers seen in the V-band light curve of ω CMa should therefore have a polarimetric counterpart. Indeed, some of the polarization variability seems to be directly related to flickering events (e.g., epochs 56,300–56,800).

Below we use the values of the field star method, because the light-curve method can have systematic errors if P_int at the epochs chosen (bracketed by the red lines in Figure 7) is nonzero. In fact, our model calculations (see Figure 8 below) indicate that this may be the case. The same can be said about the field star method, of course, as the field star may not probe the same P_IS as the target star. Our choice of the field star method is thus based solely on the fact that one method (the light-curve one) very likely suffers from systematic errors, while for the other method this is unknown.

Zoom In Zoom Out Reset image size

4.1.4. Model–Data Comparison

In this section our comparison of the model with photometric data in wavelengths other than the visible wavelength is presented.

Figure 9 shows the synthetic light curves of ω CMa. The top panel of Figure 9 illustrates the V-band modeling presented in G18. The second panel shows the synthetic UBV bands together. Interestingly, the U-band light curve shows the largest variations, which is expected as this band should be more sensitive to the inner-disk conditions, where the density varies widely (see Figure 11 in G18). The model predicts a complex behavior. In general, we see that the longer the wavelength, the slower the rate of magnitude variations in the model light curve. This is explained by the fact that larger disk volumes (from which long-wavelength continuum fluxes originate; see Figure 7 of Rivinius et al. 2013b) respond more slowly to variations in the inner disk. One interesting feature of the models is that each subsequent dissipation reaches a lower flux level when compared to the previous. This reflects the finding in G18 that in ω CMa a true quiescence value is never realized, but rather the star transitions between high and low mass-loss rate states.

Zoom In Zoom Out Reset image size

The left panel of Figure 10 shows a comparison between the observed UBV data and the color indices with the model. The model fits the data generally well. There are a few outliers, which are likely explained by the fact that we did not model the short-term flickering events (i.e., events shorter than 2 months; see G18) of the light curve as mentioned previously. A systematic mismatch, however, is seen in the color indices, most notably at B − V.

Zoom In Zoom Out Reset image size

The Strömgren LTPV data (uvby bands) are illustrated in the right panel of Figure 10. The general behavior is similar to that of the UBV data. The general shape of the curve, as well as the colors, is however quite well reproduced by the model.

The rough agreement between the VDD model and the uvby-band data is a significant result, as these bandpasses probe slightly different regions of the disk. In general the flux level is related to the disk density, while the colors probe the density gradients. Our results indicate that the density scale of the inner disk is well reproduced by the model, but the mismatch in the colors may point to inaccuracies in the density gradient. This is not surprising, given the simplistic nature of our assumptions for the disk inner-boundary conditions (see G18 and Rímulo et al. 2018 for more details).

Figure 11 shows a comparison between the observed JHKL magnitudes and colors and the model. The top panel reveals that the model is consistently brighter than the data. This discrepancy is maximum for the L band and minimum for the J band. Thus, the longer the wavelength, the larger the discrepancy. The middle and bottom panels show that the color indices are also, in general, systematically shifted with respect to the models. It should be noted that it is unlikely that the discrepancies seen for the JHKL bands are due to uncertainties in the inclination angle. For instance, changing the inclination angle to 18° in the models would cause a magnitude increment of only ≈0.02.

Zoom In Zoom Out Reset image size

Although the scarcity of IR data makes a firm conclusion difficult, it still appears that the JHKL color indices are more or less well reproduced by the VDD model, while the actual magnitudes are not. The model is always systematically brighter in the JHKL bands than the observations, which means that the model may be too dense in the outer-disk regions. This might suggest a larger α in the outer disk, which would drain it faster, making it less dense. This is discussed in Section 5.

Figure 12 displays the model and data in the VISIR Q1 (∼16.7–18.0 μm) and Q3 (∼19.1–20.0 μm) bands. Unlike the JHKL magnitudes, the data and models agree within the errors. Because the observational errors are large and only two points of the Q1 and Q3 data in total are available, these data do not allow us to confirm or disprove the tendencies seen for the JHKL bands.

Zoom In Zoom Out Reset image size

We have a comprehensive data set of spectra for ω CMa covering cycles 3 and 4 (Figure 2(c)), as well as Q2. As examples, we show a comparison between the VDD model and the observed equivalent width (EW), E/C, and peak separation (PS) for Hα and Hβ in Figure 13. The results for Hγ and Hδ are qualitatively similar and are shown in Appendix C and in Figure 19.

Zoom In Zoom Out Reset image size

These spectra have a good S/N and medium-to-high resolution. Thus, the uncertainties of the quantities are small. The uncertainties for the line-profile E/C and EW are dominated by their continuum determination, while for the PS the resolution is the limiting factor.

It is important to emphasize that we did not adjust the model to obtain an optimum fit. We used the same model and scenario described in G18 to calculate these profiles, in order to evaluate how this model performs as compared to multi-technique data.

The second and third panels of each plot in Figure 13 show the EW and E/C of the line. The EW and the E/C of a line reveal a complex interplay between the line emission and the adjacent continuum emission. The Hα line emission comes from a large volume of the disk, which responds very slowly to changes in the disk feeding rate. Conversely, the adjacent continuum responds very quickly to these changes. Therefore, when the continuum emission rises (e.g., during an outburst), the EW initially drops in magnitude and the E/C falls as well. However, when the continuum emission drops (e.g., during quiescence), the EW will increase in magnitude and the E/C will increase. These effects are more moderate for Hβ since at 4861 Å the adjacent continuum displays a much smaller range of magnitude variation along a given cycle than at 6562 Å (Figure 9).

Figure 13 can be interpreted with the above scenarios in mind. In all quiescence phases, the EW increases in magnitude (becoming more negative) and the E/C increases, as a result of the quick suppression of the inner disk, which causes the emission in the adjacent continuum to drop quickly. This initial dissipation of the inner disk does not affect the line emission. Only much later in the dissipation, when the entire disk empties, does the line emission drop. Then, the EW decreases in magnitude and the E/C also drops. At outburst, the converse happens: the inner disk fills up quickly, giving rise to a sudden increase in the continuum out to the IR. As a result, the EW decreases in magnitude and the E/C decreases (recall the apparent contradiction mentioned at the end of Section 2 and also see Figure 3).

The PS was computed by fitting a Gaussian curve to each emission peak, in order to determine its height as compared to the adjacent continuum (which was normalized to 1), as well as its velocity. The low inclination angle of ω CMa causes an almost single-peak profile in Hα whose flux comes mainly from the larger part of the disk (in comparison to the other hydrogen lines), where the Keplerian velocities are lower. Also, some of the data are of low resolution, which makes our analysis difficult. For this reason we show the PS in Figure 13 only for the spectra with a clear double-peaked structure. Comparison with the model reveals that, as is seen for the EW and E/C, there is better agreement during the outburst phases than during dissipation.

In general, the results for Hβ are similar to those for Hα. The EW curve is qualitatively reproduced, but a quantitative comparison fails mainly during the quiescence phases. Of particular significance is the close match between the data and the model for the fast decline in EW of O4. The E/C is also well reproduced. It is important to recall that since the Hβ opacity is lower than the Hα opacity, the formation volume of this line is smaller (Carciofi 2011). This can be seen by the PS values, which lie at ≈40 km s⁻¹ for Hβ, while for Hα they are, in general, smaller than 20 km s⁻¹. The larger PS indicates that Hβ is indeed formed closer to the star, where the rotational velocities are higher. The fact that the model reproduces this behavior is a significant result.

The model can reproduce these variations qualitatively, but not quantitatively. After the quick increase in magnitude of the EW at the onset of dissipation, the observed EW decreases in magnitude at a much faster rate than the model EW (the same is observed with the E/C). Since the predictions from G18 fit the visible- and IR-band light curves well, the problem likely lies in the outer disk. It appears that while G18's model predicts the correct rate of density variation in the inner disk, the corresponding rates in the outer part are too slow. In other words, the outer disk is not being drained of material fast enough. Further support to this comes from polarimetry. Recall that the observed rate of polarimetric variation is higher than the model calculations, also indicating faster emptying than in the model.

In the following, we provide some tentative explanations for this mismatch between the model and the data. One way to achieve faster dissipation rates in the outer disk is to have larger values of the viscosity; this could happen either because the temperature rises with the radius (which is not physically justified) or because the α parameter increases with distance from the star. Therefore, this might be the first hint of a radially varying α in a Be star. Another possibility is to consider that there is an unknown binary companion truncating the disk at radii smaller than the 1000 R_eq assumed here (note that this is an arbitrary assumption for the disk size in G18). If this were the case, the mass reservoir of the outer disk would be smaller, and the whole disk would dissipate faster, as suggested by the observations. Finally, a third possible explanation for the mismatch is radiative ablation. In the absence of active feeding, ablation could act in addition to viscosity to dissipate disk material. However, the results of Kee et al. (2016a and subsequent papers) indicate that ablation is more efficient at clearing the inner disk than at clearing the outer. Therefore, ablation, if included in our models, would likely make the mismatch between the rates of dissipation of the inner and outer disks worse. In the next section, we investigate the first two possibilities, namely larger α in the outer disk and binary truncation. Also, we discuss the influence of ablation in more detail.