The Be Star 66 Ophiuchi: 60 Years of Disk Evolution

Classical B-emission (Be) stars are rapidly rotating, main-sequence, B-type stars that form outwardly diffusing disks of gas that have been ejected from the stellar surface. These geometrically thin disks are composed of almost fully ionized hydrogen gas and rotate around the stellar equator in Keplerian fashion (Rivinius et al. 2013). The rapid rotation of these stars (Slettebak 1982) is believed to lead to the disk formation, which is likely assisted by nonradial pulsations (Rivinius et al. 2003; Baade et al. 2016). The disks are characterized by infrared and radio excesses in the stellar spectra, polarized light, and the presence of hydrogen series emission lines (Rivinius et al. 2013).

Many Be stars show variability over timescales ranging from minutes (Goraya 2008) to decades (Okazaki 1997). The rapid variability observed from these systems is commonly associated with localized mass ejections (Balona 1990), β-Cephei type pulsations (Balona & Rozowsky 1991), and nonradial pulsations (Baade 1982; Rivinius et al. 2003). Some disks exhibit cycles of growth and dissipation, which are dependent on the viscosity of the gas as outlined in the viscous decretion disk (VDD) model (Lee et al. 1991; Papaloizou et al. 1992; Klement et al. 2015). This type of variability occurs typically over periods of decades (e.g., Wisniewski et al. 2010).

The star 66 Ophiuchi (HR 6712, HD 164284) is a bright (m_v ∼ 4.8 mag), multiple star system (Štefl et al. 2004) containing a classical Be star of spectral type B2Ve (Floquet et al. 2002), at a distance of ∼199.6 pc. ³ Since 1957, observations indicate that 66 Oph built a large disk that it subsequently lost over a period of dissipation that started around 1990 and finished about 20 years later. These events make 66 Oph an ideal system for studying the physics of disk dissipation.

Many observational campaigns have focused on this star because of the variability of its disk. In 1957, the transition of the Hα line into emission signaled the onset of disk formation (Rakotoarijimy & Herman 1958). Grady et al. (1987), Percy & Attard (1992), and Percy & Bakos (2001) found that metal lines at UV wavelengths showed evidence of recurrent episodic mass loss from 1982 to 1985, and again from 1989 to 1999, using UBV photometry. Hanuschik et al. (1996) reported that the ratio of the violet (V) and red (R) peaks of the Hα profile began to vary in 1988, and continued through 1995 when it had a period of five years. Štefl et al. (2004) showed that the V/R ratio was constant (i.e., V = R) at the onset of dissipation, suggesting the disk was axisymmetric. Floquet et al. (2002) gave a detailed history of the pulsation period of 66 Oph before the period of disk dissipation, which they claimed began in 1992. Since then, there has been a slow decline of the Hα emission strength with a transition to an absorption line in 2010 (Sabogal et al. 2017). The line profile has remained unchanged since then.

In this study, we investigated the evolution of 66 Oph's disk at the onset and throughout the dissipation period. We characterized the physical state of the disk before dissipation and the dynamics of the disk during dissipation. This includes identifying the density structure and radial extent of the disk, as well as the viscosity and temperature distribution. This was accomplished by finding the best model to reproduce observations of the spectral energy distribution (SED), Hα spectral line, V-band polarization, and V-band photometry during dissipation. By including observations from previous works in addition to presenting new observations, our unique data set allowed us to model the complete dissipation of 66 Oph's disk for the first time.

The models presented in this work were created following similar methods to Haubois et al. (2012), Carciofi et al. (2012), Rímulo et al. (2018), and Ghoreyshi et al. (2018). The 1D dynamical code SINGLEBE was used to compute the disk surface density of an axisymmetric disk as a function of time, given a disk viscosity. The resulting disk density distributions were then input to the 3D radiative transfer code HDUST to calculate the emergent spectrum. In Section 2, we describe the observations compiled over the period of 1957 to 2020. Section 3 gives a detailed description of our method to determine the parameters of the star. Section 4 describes our modeling routines along with the results of the modeling. Section 5 provides a discussion and summary of our work including a comparison to other findings in the literature.

We compiled observations from a variety of different sources to investigate the disk's evolution over many decades. Figure 1 shows the V-magnitude photometry, Hα equivalent width (henceforth, EW), and continuum V-band polarization from the onset of disk growth and subsequent dissipation from 1957 to 2020. We adopt the convention that negative EW means flux above the continuum.

Zoom In Zoom Out Reset image size

Most of the observations used in this work were from the literature. We used V-magnitude photometry from the works of Haupt & Schroll (1974), Pavlovski et al. (1997), Percy et al. (1997), and Floquet et al. (2002). We also acquired Hα EWs from the works of Slettebak & Reynolds (1978), Ghosh (1988), Krishnamurthy (1999), Hanuschik et al. (1995), Floquet et al. (2002), and from the Be Star Spectra Database (BeSS). ⁴ Reported values of the Hα EW and V-band polarization were acquired from the archive for the Lyot Spectropolarimeter ⁵ and the Halfwave Spectropolarimeter (HPOL) at the University of Wisconsin–Madison Pine Bluff Observatory, which were reduced in the work of Draper et al. (2014). We also used observations of the V-band polarization published by Poeckert et al. (1979) and Hayes (1983).

A number of previously unreported observations of 66 Oph are also presented. We used measurements of the Hα EW determined from observations of the Hα spectra made using the fiber-fed échelle spectrograph attached to the 1.1 m John S. Hall telescope at the Lowell Observatory in Flagstaff, Arizona. Observations from this instrument were obtained between 2005 and 2007, with a resolving power of R = 10,000. The reduction process of these observations follows that previously described in Jones et al. (2017).

Our models are constrained by observations of the Hα EW and V-band polarization, acquired respectively using the MUSICOS spectrograph and the IAGPOL polarimeter at the Pico dos Dias Observatory (OPD), operated by the National Astrophysical Laboratory of Brazil (LNA) in Minas Gerais, Brazil. These observations were reduced with packages developed by the Beacon group ⁶ and described in Magalhaes et al. (1984, 1996) and Carciofi et al. (2007).

The most recent Hα EW observation was obtained by the NRES spectrograph at the Las Cumbres Observatory LCOGT network. Details about the instrument and reduction process of this observation can be found in Brown et al. (2013). Since this observation was acquired while 66 Oph was diskless, it was used for comparison to our diskless models. This observation is later presented in Section 4.2.

Over a 63 yr period, the V-band photometry, m_ν , was observed to range between 4.5 mag < m_ν < 4.9 mag. During dissipation, the episodic variability continued with nine outbursts of between 1991 and 2008, while the overall brightness asymptotically dimmed.

As the V-band photometry available to us is sparse during the disk building phase, the variation of the light could not be used to model the evolution of the disk as was done in the studies of Ghoreyshi et al. (2018) and Rímulo et al. (2018). Photometric observations in other bands (UV, IR, etc.) are also sparse, with only snapshots available. Here, we use the V-band photometry along with observations of the entire SED, the Hα line profile, and the percentage of polarized light over time.

The vertical, solid gray line in Figure 1 indicates the onset of dissipation (and is further discussed in Section 4.1). As the dissipation event begins, the continuum flux drops and the inner disk re-accretes causing the Hα EW to increase. This is seen in Figure 1 around 1995, as shown by the observations from Krishnamurthy (1999) and Hanuschik et al. (1995). After this time, the EW began to steadily decrease until the line went into absorption in 2010 (indicated by the vertical, dashed gray line). Our most recent observation, obtained by NRES/LCOGT in 2020 March, indicates that 66 Oph continues to be diskless.

The bottom panel of Figure 1 shows the change in the observed polarization, p_v , with time for 66 Oph. Since 1980, the percent polarization slowly decayed until it approached the interstellar base level. From the base level polarimetric observations obtained by OPD/LNA in 2017 (listed in Table 1), the V-band interstellar polarization is ∼0.63% with a polarization position angle of θ ≈ 857. This average value was subtracted from each of the observations using the scheme outlined by Quirrenbach et al. (1997). The intrinsic polarization of 66 Oph obtained is discussed further in Section 4.2.2.

The star's SED was also used to refine our models. We selected 112 observations of the UV portion of 66 Oph's SED observed by the International Ultraviolet Explorer (IUE) telescope from the INES database (González-Riestra et al. 2001) following the selection procedure described by Freire Ferrero et al. (2012). The data selected were observed using the large aperture and high dispersion modes on the LWR and SWP cameras to ensure proper flux calibration and high spectral resolution. The observations were obtained over the period of 1980 March to 1995 September, during which the observed flux showed no significant changes. We chose to remove the IUE data beyond 0.3 μm due to instrumental limitations that cause significant uncertainty (González-Riestra et al. 2001).

Visible and IR wavelength observations were obtained from the CDS Portal application from the Université de Strasbourg. Observations at these wavelengths were acquired while the disk was present (e.g., IRAS in 1989) and absent (e.g., AKARI in 2006–07 and WISE in 2010), providing upper and lower limits for modeling the SED. In addition to indicating the onset of dissipation and transition to absorption, the gray lines in Figure 1 also lie at the approximate dates of observation of IRAS (solid gray line) and WISE (dashed gray line) data. We also collected radio observations from the Very Large Array (VLA) and ATCA (Apparao et al. 1990; Clark et al. 1998) and JVLA (Klement et al. 2019). These observations are presented alongside our models in Section 4.1.

2.1. The Interstellar Polarization

Figure 2 shows a Stokes QU diagram of the V-band polarization obtained by HPOL and by OPD/LNA. The observed path on the QU diagram during the disk dissipation phase is upward and to the right. As discussed by Draper et al. (2014) and more recently by Ghoreyshi et al. (2021), the process of formation and/or dissipation of a Be star disk is associated with a linear trend in the QU diagram, during which the polarization level rises or lowers in response to changes in the disk density, but the polarization angle remains steady. By fitting a linear regression to the HPOL observations, and fixing it to the assumed diskless 2016 observation, we determined the polarization position angle of the disk to be ∼98° ± 3°. This was also confirmed by fitting a Serkowski law (Serkowski 1973) to the OPD/LNA observations in Table 1, using

$\begin{eqnarray}&&P(\lambda )/P({\lambda }_{\max })=\exp [-1.15\,{\mathrm{ln}}^{2}({\lambda }_{\max }/\lambda )].\end{eqnarray} \tag{ 1 }$

Zoom In Zoom Out Reset image size

We determined this law to best fit the observations in Table 1 when $P({\lambda }_{\max })$ = 0.62 ± 0.02% at 0.61 ± 0.04 μm. Subtracting this spectrum from the HPOL observation, we accurately modeled the polarization of the disk (shown later in Section 4.2.2). From this spectrum, the wavelength averaged polarization position angle of the disk was computed to be 96° ± 4°, in perfect agreement with the estimate made using the QU diagram.

The use of Markov Chain Monte Carlo (MCMC) routines along with Bayesian statistics has recently found success in the modeling of Be stars (e.g., Rímulo et al. 2018; Mota 2019; Suffak et al. 2020; B. C. Mota et al. 2021, in preparation). In this Section, we use these methods to determine which stellar parameters best reproduce observations from the IUE, previously described in Section 2. The UV spectrum is not strongly affected by the disk, so we used it for our fitting procedure to find the stellar parameters. We determine values for the stellar mass M, critical fraction of rotation W (as defined in Equation (6) of Rivinius et al. 2013), age t/t_ms (where t_ms is the main-sequence lifetime), inclination i, distance d, and the degree of interstellar reddening E(B − V), and from these compute the derived parameters listed in Table 2.

Table 2. The Best-fitting Stellar Parameters for 66 Oph Computed with emcee, and Derived from the Computed Values Using the Models of Georgy et al. (2013)

Best-fit		Derived
Parameters	Values	Parameters	Values
M [M⊙]	${11.03}_{-0.53}^{+0.55}$	L [L⊙]	${8200}_{-1300}^{+1600}$
W	${0.52}_{-0.05}^{+0.05}$	Teff [K]	${25940}_{-750}^{+800}$
t/tms	${0.33}_{-0.08}^{+0.10}$	$\mathrm{log}g$	${4.17}_{-0.06}^{+0.05}$
i [°]	${57.6}_{-6.8}^{+6.8}$	Rpole [R⊙]	${4.50}_{-0.23}^{+0.32}$
d [pc]	${208.6}_{-9.0}^{+8.9}$	Req [R⊙]	${5.11}_{-0.26}^{+0.37}$
E(B − V)	${0.22}_{-0.01}^{+0.01}$	$v\sin (i)$	${290}_{-9}^{+11}$

Download table as:  ASCIITypeset image

We use a grid of 770 diskless Be star models, called BeAtlas, to create a parameter space for evaluating the observed spectrum from IUE. The BeAtlas grid was originally computed by Mota (2019) using a 3D non-LTE Monte Carlo radiative transfer code called HDUST (Carciofi & Bjorkman 2006). This code simulates Be stars with or without disks by solving the coupled problems of radiative equilibrium, statistical equilibrium, and radiative transfer to compute synthetic observables. Table 3 summarizes the ranges of the stellar parameters for the BeAtlas models; the masses are consistent with those computed by Georgy et al. (2013).

The parameter space was sampled using emcee (Foreman-Mackey et al. 2013), a Python code of the Affine Invariant MCMC Ensemble sampler (Goodman & Weare 2010). We defined Gaussian prior distributions with mean and variance taken from literature values (see Table 4). We used the parallax from Hipparcos (van Leeuwen 2007), since the parallaxes from GAIA's DR2 catalog (Gaia Collaboration et al. 2018) contained large errors for bright stars, including 66 Oph. Recently, we found recomputing the stellar parameters using the parallax reported in GAIA's eDR3 catalog (Gaia Collaboration et al. 2020; ∼4.90 ± 0.37 mas, or ∼${204}_{-14}^{+17}\,\mathrm{pc}$) produced a consistent set of parameters.

Before fitting, each model was also artificially reddened using the standard Fitzpatrick (1999) parameterization, with E(B − V) as a free parameter in determining each model's goodness of fit. We used a $\mathrm{log}({\chi }^{2})$ likelihood function to evaluate the fit of the models to the observations.

Our emcee computation used 30 walkers and 50,000 steps with a burn-in of 5000 steps, which were chosen by following the guidelines discussed by Foreman-Mackey et al. (2013). Approximately 27% of the steps proposed by the MCMC were accepted, which is consistent with the range recommended by Foreman-Mackey et al. (2013), when using a walker step size of 3.0.

Figure 3 shows the resulting probability density functions (PDFs) for the stellar parameters and correlation between these parameters as predicted by emcee. The top-right corner of Figure 3 shows UV data along with a model computed with the best-fit parameters in red. The gray lines show the model predictions corresponding to the last step of each walker. The residuals between the data and the model are shown directly below this panel.

Zoom In Zoom Out Reset image size

The stellar parameters predicted by emcee are used to compute the mean stellar luminosity L, mean effective temperature T_eff, mean surface gravity g, polar radius R_pole, equatorial radius R_eq, and $v\sin (i)$ by interpolating the Geneva stellar evolutionary models from Georgy et al. (2013; more details are available in Mota 2019). Table 2 lists the most probable stellar parameters and the additional parameters derived from the Geneva models. The mass, temperature, and radius are consistent with standard values for B0 to B2 stars (Cox 2000), and are in agreement with Waters et al. (1987), Floquet et al. (2002), Frémat et al. (2006), and Vieira et al. (2017). Values of $v\sin (i)$ reported in the literature range from 200 km s⁻¹ (Bernacca & Perinotto 1970) to 292 km s⁻¹ (Floquet et al. 2002); our derived value agrees more closely with the upper limit. Our stellar parameters are further confirmed in Section 4.1, by comparison with the observed diskless SED.

In Section 4.1 we determine the quasi-steady state of the disk by SED fitting and Hα line modeling prior to the dissipation. We model the Hα observation obtained by the HPOL Be star campaign on 1989 June 28. This observation (previously discussed by Draper et al. 2014) corresponds to our determination of the onset of disk dissipation in 1989 March. This observation was obtained at a key time when 66 Oph's disk was neither growing nor dissipating. The evolution of the disk before and after the quasi-steady state is modeled in Section 4.2. These models consider four different scenarios: whether the disk is isothermal or nonisothermal, and whether the disk builds to a steady state then dissipates, or dissipates from the steady state. These scenarios are presented as follows: Section 4.2.1—isothermal disks that build then dissipate, Section 4.2.2—isothermal disks that dissipate only, Section 4.2.3—nonisothermal disks that build then dissipate, and Section 4.2.4—nonisothermal disks that dissipate only. Models for each scenario were constrained by SED and Hα line fitting, and further compared to polarization when appropriate.

4.1. Modeling the Quasi-steady State Phase

In the VDD model, gas is injected into the inner disk and diffuses outward via viscous forces (Porter 1999). VDDs were initially modeled by Porter (1999) and later by Okazaki (2001), Bjorkman & Carciofi (2005), Jones et al. (2008), and Lee (2013). These studies showed that with a constant mass injection rate the disk's density profile becomes relatively constant as a steady state is reached over months, years, or longer, since gas that is moving outward is replenished at the innermost boundary. However, Haubois et al. (2012) claims that a true steady state is never realized but that a quasi-steady state phase is observed in many Be stars (e.g., Klement et al. 2015).

The Hα EW did not vary significantly between 1985 and 1995 (as seen in the middle panel in Figure 1). During this time, 66 Oph's disk was approximately quasi-steady; however, the drop in V-band magnitude and polarization, which form closer to the star (Faes et al. 2013), suggests that the dissipation began between 1989 and 1990 (indicated by the vertical, solid gray line in Figure 1).

Our initial goal is to quantify the density structure of 66 Oph's disk while in this quasi-steady state so that we can use this model as a baseline to follow subsequent dissipation.

Often in the literature, a power-law density distribution is used to represent a snapshot of the disk in time. We use the density structure originally adopted by Waters (1986) for IR continuum observations,

$\begin{eqnarray}&&\rho (r)={\rho }_{0}{\left(\displaystyle \frac{{R}_{\mathrm{eq}}}{r}\right)}^{n}.\end{eqnarray} \tag{ 2 }$

Here, ρ₀ is the density at the innermost disk radius in the equatorial plane, r is the radial distance from the central star, and n is the power-law exponent.

We computed a grid of 61,568 models with densities from ρ₀ = 10⁻¹⁰ g cm⁻³ to 10⁻¹³ g cm⁻³ for every quarter magnitude, with n = 2.0 to 3.5 in steps of 0.1, inclination from i = 0° to 90° in steps of 25, and outer disk radii of R_out = 20, 30, 40, 50, 60, 80, 100, and 200 R_eq. For each model, HDUST (previously described in Section 3) was used to compute the SED and Hα line profile.

Noncoherent electron scattering within the disk also affects the shape of line profiles by broadening the wings of the line (Hummel & Dachs 1992). This process is not accounted for in HDUST, so to approximate the effect, we assume that a fraction f of the Hα flux is scattered by electrons with thermal velocities ranging from v_e = 300 km s⁻¹ to 800 km s⁻¹, and that the scattered flux has a Gaussian profile. Therefore, a new line profile F_new is obtained by

$\begin{eqnarray}&&{F}_{\mathrm{new}}(\lambda )=(1-f){F}_{\mathrm{nc}}\lambda +f\times {F}_{{\rm{c}}}(\lambda ),\end{eqnarray} \tag{ 3 }$

where F_nc is the nonconvolved line profile predicted by HDUST and F_c is the convolution of F_nc with a Gaussian with FWHM of v_e . The convolved line profiles were then fit to the HPOL 1989 spectrum using emcee (unique from the routine which used emcee previously described in Section 3). The fits were evaluated by interpolating the model fluxes to the wavelengths of the observed data and then calculating ${\chi }_{\nu }^{2}$ values.

After convolving, we found a subset of models within 1σ of the observed Hα EW (−58 Å). The best-fit model according to ${\chi }_{\nu }^{2}$ testing has a density of ρ₀ = 2.5 × 10⁻¹¹ g cm⁻³, with n = 2.6, at i = 575, and R_out = 100 R_eq, with a fit of ${\chi }_{\nu }^{2}=1.17$ when using F_c = 0.22 and v_e = 700 km s⁻¹. Figure 4 shows the observed profile (black line), and the best-fit model before (dotted line) and after (dashed line) convolving. These profiles all have EW's of −58 Å. The thermal velocity used for the convolution corresponds to a kinetic temperature of T_K = 10,100 K. This agrees with the disk temperature of ∼10,500 K at the edge of the optically thick Hα emitting region at ∼10.1 R_eq (computed by following Appendix D of Vieira et al. 2017). This process of convolving lines was previously used by Klement et al. (2015), without MCMC fitting, for the Be star β CMi. β CMi is a late-type star with a spectral type of B8Ve and is expected to have a cooler disk. We note that these authors obtained a value of f = 0.6 and a smaller value of v_e = 300 km s⁻¹, consistent with its disk temperature.

Zoom In Zoom Out Reset image size

In the case of a disk with n = 3.5, the Hα emitting volume probes a significant portion of the disk (Tycner et al. 2005) typically thought to be the innermost 20 R_eq (Faes et al. 2013). However, we found these regions become increasingly extended as n decreases. This means that for an early-type star like 66 Oph, lower values of n result in more ionized material at larger distances. Figure 5 shows the Hα EW and corresponding ${\chi }_{\nu }^{2}$ value of the best-fit model with increasing R_out. Here, with n = 2.6, we find the EW increases until the disk reaches a radius of ∼100 R_eq, beyond which there is no significant increase. Similarly, ${\chi }_{\nu }^{2}$ improves as R_out increases and the Hα flux increases until ∼100 R_eq, with no appreciable change at larger radii. This provides an upper limit on the size of the Hα emitting region.

Zoom In Zoom Out Reset image size

The star's SED (Figure 6) was also used to independently constrain our model values of ρ₀ and n. The SED upper bound, prior to dissipation, is set by IRAS (1989) in the IR, as well as radio observations from VLA (1988). Conveniently, the upper bound from IRAS was observed just before the onset of dissipation, providing important observational constraints for the quasi-steady state phase. Data from ATCA (1997) and JVLA (2010) are also on the upper bound, as the disk dissipation had not significantly changed the disk's radio emission, which comes from the outermost disk. The SED's lower bound is constrained by observations from WISE (2010) and AKARI (2006–07). Observations of the lower bound were obtained after dissipation and likely correspond to a diskless state. Using the diskless model of Section 3, we find a very reasonable fit to the observations at visible to radio wavelengths observed after 1989 with ${\chi }_{\nu }^{2}=1.12$.

Zoom In Zoom Out Reset image size

We fit the upper-bound SED with the same grid of models used above for Hα fitting. We find that a model of ρ₀ = 10⁻¹¹ g cm⁻³, n = 2.4, i = 575, and R_out = 100 R_eq (shown as the solid black line in Figure 6, and the black dot in Figure 7) produced the best fit to the visible, IR, and radio photometry observed before 1989 with ${\chi }_{\nu }^{2}=1.09$. Increasing the outer disk radius to 200 R_eq marginally improved the fit to ${\chi }_{\nu }^{2}=1.07$. The model that best fit the Hα line profile (ρ₀ = 2.5 × 10⁻¹¹ g cm⁻³ and n = 2.6) fit the SED observations before 1989 with ${\chi }_{\nu }^{2}=1.15$ for R_out = 100 R_eq, or ${\chi }_{\nu }^{2}=1.14$ for R_out = 200 R_eq.

Zoom In Zoom Out Reset image size

Figure 7 shows a comparison of the best-fit models to the Hα line profile and SED according to the ${\chi }_{\nu }^{2}$ fitting. The contours show the models that fit best to within 1σ. The overlap has a range of ρ₀ = 8.5 × 10⁻¹² g cm⁻³ to 3 × 10⁻¹¹ g cm⁻³, and n = 2.4 to 2.7. Within the same overlapping region, only the ρ₀ = 2.5 × 10⁻¹¹ g cm⁻³ and n = 2.6 model satisfies both the SED and Hα line profile fitting. The model that best fit the SED produced a relatively poor fit to the Hα line profile, with ${\chi }_{\nu }^{2}=3.2$.

The ρ₀ = 2.5 × 10⁻¹¹ g cm⁻³ and n = 2.6 model predicts a V-band magnitude of 4.72 mag, which is consistent with the range of 4.76–4.55 mag observed in 1989 (Floquet et al. 2002). The same best model predicts the V-band polarization to be ∼0.72%, while observations obtained 5 years prior to dissipation show the V-band polarization to range from 0.58% to 0.65% (after subtracting the interstellar polarization, as described in Section 2).

Since the ρ₀ = 2.5 × 10⁻¹¹ g cm⁻³ and n = 2.6 model with i = 575 and R_out = 100 R_eq is clearly successful in reproducing the observed Hα line profile, SED, V-band magnitude, and V-band polarization, we conclude it is the best representation of the quasi-steady state phase for 66 Oph's disk.

4.2. Hydrodynamic Calculations

In this section, we model the formation of 66 Oph's disk and the subsequent 21 years of dissipation until the Hα emission went into absorption. We use the best-fit quasi-steady model from Section 4.1 as a midpoint between the building and dissipation. To simulate the disk evolution, we use the 1D dynamical code SINGLEBE, developed by Okazaki et al. (2002) and later used by Carciofi et al. (2012), Haubois et al. (2012), Rímulo et al. (2018), and Ghoreyshi et al. (2018).

SINGLEBE solves the 1D fluid hydrodynamic equations (Pringle 1981) for the surface density of a viscous disk using the thin disk approximation. The disk is assumed to be axisymmetric, Keplerian, and in vertical hydrostatic equilibrium. See Okazaki et al. (2002) and Rímulo et al. (2018) for more details on SINGLEBE.

The evolution of VDDs is controlled by kinematic viscous forces within the disk. SINGLEBE adopts the commonly used α-prescription for viscosity, defined by Shakura & Sunyaev (1973) as

$\begin{eqnarray}&&\nu (r)=\displaystyle \frac{2}{3}\alpha (r){c}_{s}(r)H(r),\end{eqnarray} \tag{ 4 }$

where ν is the disk viscosity, r is the radius in the disk, and H is the disk scale height. The sound speed is given by

$\begin{eqnarray}&&{c}_{s}(r)={\left[({k}_{B}T(r))/(\mu {m}_{a})\right]}^{1/2},\end{eqnarray} \tag{ 5 }$

where k_B is the Boltzmann constant, T is the gas kinetic temperature, μ is the mean molecular weight, and m_a is the atomic mass unit. Note, for a nonisothermal disk, the radial dependence of the temperature can change the sound speed, which in turn changes the viscosity in Equation (4). In this case, α(r) and T(r) become degenerate and indistinguishable. If one parameter were to be held constant, the other parameter may be adjusted to produce the same viscosity. Therefore, we distinguish isothermal disks from nonisothermal disks based on constant or variable α T with radius, respectively.

Viscosity affects disk evolution by influencing the outflow rate of the gas (Lee et al. 1991). Given α(r), a target inner surface density Σ(r = R_eq); (i.e., Σ₀) or inner volume density ρ₀, the outer disk radius R_out, and a value for the mass injection rate at the base of the disk ${\dot{M}}_{\mathrm{inj}}$, SINGLEBE computes the surface density as a function of time and distance from the star, Σ(r, t). Details about the computational procedure and assumed boundary conditions can be found in Rímulo et al. (2018).

Using Equation 3.6 of Okazaki (1991), the mass density can be computed as

$\begin{eqnarray}&&\rho (r,t)=\displaystyle \frac{{\rm{\Sigma }}(r,t)}{\sqrt{2\pi }H(r,t)}.\end{eqnarray} \tag{ 6 }$

We note that, assuming a constant α T, hydrodynamic theory for an axisymmetric disk with Keplerian rotation in a quasi-steady state predicts that the mass density will have approximately a radial power law (Equation (2)) with n = 3.5 (Bjorkman & Carciofi 2005).

We input the stellar parameters determined in Section 3 (see Table 2), ρ₀, n, and R_out of the best-fit model from Section 4.1, and the duration of the growth and dissipation events to SINGLEBE. We use α = 1.0 to ensure rapid evolution (Haubois et al. 2012); however, as α can be scaled by changing the length of the evolutionary epoch (as illustrated in Figure 1 of Haubois et al. 2012), our models can explore all values of α. For each model, the disk's innermost radius was set to 1 R_eq and the radial grid points were spaced logarithmically from the inner radius. Following Rímulo et al. (2018), we define the onset of dissipation to occur at 0 yr for all models. The density profiles from SINGLEBE were input to HDUST to compute the SED, Hα line profile, and polarization at times corresponding with the observed data.

4.2.1. Building and Dissipating an Isothermal Disk

Here, we investigate an isothermal disk that builds to a steady state and then dissipates, which we refer to as the constant α T, n = 3.5 scenario. The isothermal assumption also makes our models similar to the mixed models used in Carciofi & Bjorkman (2008), where the disk is isothermal in the fluid equations and variable in temperature when computing observables using HDUST.

To match the base density found in Section 4.1, this scenario requires a mass injection rate of $\dot{M}=8.5\times {10}^{18}\,{\rm{g}}{{\rm{s}}}^{-1}$ (or ∼1.3 × 10⁻⁷ M_☉ yr⁻¹) at R = 1.02 R_eq, of which ∼99% falls back onto the star. We chose a constant disk temperature at 60% of T_eff (from Table 2), following Carciofi & Bjorkman (2006).

Figure 8(a) shows density profiles of this scenario after the building period during the dissipation. The quasi-steady state is indicated by the thick, black line on the panel. We note that it has the expected value of n = 3.5, which differs from the n = 2.6 of the best-fitting quasi-steady model. During dissipation, the disk empties in an "inside-out" fashion, which is expected for Be stars, as discussed by Poeckert et al. (1982), Rivinius et al. (2001), and Clark et al. (2003). After the dissipation, the maximum density of the disk is ∼5 × 10⁻¹³ g cm⁻³.

Zoom In Zoom Out Reset image size

Figure 9 shows the Hα line profiles, the disk temperature profiles, and the visible, IR, and radio SEDs compared to observed data (previously shown in Figure 6), at different times during dissipation. The SEDs in the three right panels show that this scenario is unable to reproduce the observed flux before dissipation. This is most notable at IR wavelengths where the flux is ∼1–2 orders of magnitude fainter. In the top-left panel of Figure 9, the Hα line at 0 yr shows a peak flux of ∼1.5 times the continuum flux. In contrast, the HPOL 1989 observation (previously shown in Figure 4) shows a peak flux of ∼11 times the continuum flux.

Zoom In Zoom Out Reset image size

The n = 3.5 disk also has a local maximum of the peak Hα flux at 3 yr after the dissipation started. This is due to an increase of stellar ionizing radiation (and corresponding increase in Hα flux) reaching larger radial distances of the inner disk due to the reaccretion of material closest to the star. (See the bottom left panel of Figure 9.) After this, the disk continues to dissipate and the Hα flux drops smoothly. This phenomenon was also observed for disk models with α = 0.1 and α = 0.5, with the flux increase occurring at later times in the dissipation.

In Figure 10 we show the SED for this scenario after 21 years of dissipation (black line), and also the best quasi-steady-state model (gray line). As the disk dissipates, the model approaches the theoretical diskless SED. This model closely fits the UV and IR observations; however, it produces a weaker radio excess than the 2010 JVLA observations.

Zoom In Zoom Out Reset image size

Figure 11 shows the Hα line profile (black line) for this scenario after 21 years of dissipation. The line profile shows no emission and is consistent with the absorption profile of the NRES/LCOGT 2020 observation. Note that the NRES/LCOGT observation closely matches with the observations from BeSS obtained in 2010.

Zoom In Zoom Out Reset image size

Figure 12 shows the Hα EW during dissipation for this scenario (black line). Here it is apparent that the disk could only build to approximately half of the observed EW of the quasi-steady state and dissipates at a much faster rate than the observations.

Zoom In Zoom Out Reset image size

In conclusion, the constant α T, n = 3.5 scenario fails to reproduce the two main observables chosen to constrain the models, namely the SED and Hα line emission. This is further illustrated in Figure 8(a), which shows that this model produces a disk far less massive than what is required. This will be further discussed below.

4.2.2. Dissipating an Isothermal Disk from a Defined Density

Here, we use an ad hoc scenario that starts to dissipate from the quasi-steady state with n = 2.6 in order to bypass the building phase. This constant α T, n = 2.6 scenario uses the same α and T_disk as the constant α T, n = 3.5 scenario.

Figure 8(b) shows that as dissipation begins, this scenario quickly moves toward a steeper density slope. During the first five years, the gas from the outer disk moves to the inner disk, preventing the density of the disk between ∼5 R_eq and 10 R_eq from depleting. As dissipation continues, the entire disk appears to dissipate at a steady rate across all radii.

Figure 10 shows the SED for this scenario (blue dashed line) closely matches the visible and IR observations after 21 years of dissipation. An excess of radio flux is maintained as gas remains in the outer disk.

Figure 11 shows no Hα emission (blue dashed line) remains after dissipation, consistent with the diskless NRES/LCOGT 2020 observation. In Figure 12, the Hα EW shows that this scenario (blue dashed line) requires ∼17 years to reproduce the diskless profile when α = 1.

Since α and R_out both affect the rate of disk dissipation, they must be determined simultaneously. Figure 13 shows an interpolated grid of models for this scenario computed over the ranges of 0.2 < α < 1.5 and 100 R_eq < R_out < 1000 R_eq. The black line indicates which models dissipate from the initial state to having no Hα emission in 21 yr. The combination of α = 0.4 and R_out = 115 R_eq was determined to best reproduce the observations acquired after 2, 5, 8, 9, 12, 14, 17, 19, and 21 years of dissipation. Figure 12 shows that this model closely reproduces the Hα EW curve. This model is further compared to the data in Figures 14–16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each stage of dissipation in Figure 14, noncoherent electron scattering was accounted for using the method described in Section 4.1. As previously mentioned, the 0 yr model required the electron thermal velocity in the disk to be 500 km s⁻¹ when 22% of the light is scattered by the electrons. The electron thermal velocity required to broaden the line at all times during the dissipation is ∼560 ± 10 km s⁻¹ (which corresponds to a kinetic temperature of ∼12,400 ± 500 K). The fraction, f, of scattered light drops during dissipation, going from f = 20% at 5 yr to f = 14% at 12 yr, and f = 4% at 19 yr. This is expected, as f should decrease with decreasing disk density.

Figure 15 shows the observed polarization from HPOL near the onset of dissipation. The constant α T, n = 2.6 scenario with α = 0.4 and R_out = 115 R_eq is overplotted for comparison at the same time intervals as in Figure 14.

Zoom In Zoom Out Reset image size

In Figure 16, we compare to each of the observables for α = 0.1, 0.4, and 1.0. The errors for each model were computed using a 1σ deviation of ten unique simulations run for each model. The EW begins to change most rapidly nine years after the onset of dissipation, and decreases steadily for the following ten years until the line drops into absorption. For both the V-band magnitude and the continuum polarization the most rapid change occurs between zero and five years. During this time, the α = 0.4 disk dims by ∼0.1 mag. The percent polarization intrinsic to the α = 0.4 disk drops steadily by ∼0.2% every five years. For comparison, the bottom panel shows the observed polarization position angle intrinsic to the disk, computed from the percent polarization. We find polarization position angles range between ∼94° and ∼104° and appear to remain constant through the dissipation event. These polarization position angles are consistent with the angle of 98° determined in Section 2.

Zoom In Zoom Out Reset image size

4.2.3. Building and Dissipating a Nonisothermal Disk

Here, we investigate both the building and dissipation phases using the assumption of a nonisothermal disk. For reasons made apparent below, we refer to these models as the variable α T, n = 2.7 scenario. In this case, the fluid equations no longer limit the density profile to n = 3.5, however we require knowledge of the now inseparable α T parameter. We assume a power-law relation between the product α T and the disk radius given by,

$\begin{eqnarray}&&\alpha (r)T(r)={\alpha }_{0}{T}_{0}{\left(\displaystyle \frac{{R}_{\mathrm{eq}}}{r}\right)}^{C},\end{eqnarray} \tag{ 7 }$

where α(r) and T(r) are the disk viscosity and disk temperature at radius r, α₀ and T₀ are the disk viscosity and disk temperature at r = R_eq, and C is a constant power-law index. We chose α₀ = 1 and T₀ = 60% of T_eff.

We tested values of C ranging from 0.1 to 100 and found C = 2 to most closely match the n = 2.6 density slope. As shown in panel (c) of Figure 8, the disk was unable to build to n = 2.6, with the majority of the outer disk having n ≈ 2.7. The timescales required to build the outer disk to n = 2.6 were on the order of 1000 yr, significantly greater than time period that 66 Oph built its disk. The long viscous timescales at the outer disk are explained by the low viscosity inferred from Equation (7). For instance, assuming a constant temperature and C = 2, α(r = 100 R_eq)/α₀ = 0.0001.

Figure 10 shows that the SED for this scenario (yellow line) reproduces the observed visible spectrum and radio excess after 21 years of dissipation; however, it remains too bright at IR wavelengths. The Hα line profile is shown in Figure 11 (yellow line) after accounting for noncoherent electron scattering using the process outlined in Section 4.1. We see that the emission remains too strong, evidence that the disk depletes too slowly, which is also seen in the Hα EW in Figure 12 (yellow line).

4.2.4. Dissipating a Nonisothermal Disk from a Defined Density

For our final dynamical scenario, we again adopt an ad hoc solution starting from the best quasi-steady model of Section 4.2.1 to study the dissipation phase using the nonisothermal assumption. We refer to these models as the variable α T, n = 2.6 scenario. In this case, we assume that the outer disk could have built through past (undocumented) disk events, since Be stars remain on the main sequence for most of their lifetime (Georgy et al. 2013). We chose to use the α T profile described in Equation (7), with C = 2.

Figure 8(d) shows that the disk clears from the "inside out." The most rapid change to the disk's density occurs within the first five years, beyond which the density decreases more slowly. The final five years of dissipation produces a change in density at 10 R_eq of ∼0.5%. As before, we find the viscous timescales of the outer disk are too large since α T decreases too rapidly with radius.

Figure 10 shows this scenario (red line) does not match the diskless SED observations as closely as the isothermal models. After dissipation, this scenario more closely matches the diskless model from the UV to the near-IR, and more closely matches the quasi-steady state in the far-IR. Similarly to the variable α T, n = 2.7 scenario, the excess radio flux remains bright, matching the observations.

Figure 11 shows Hα emission (red line) remains too strong after 21 years of dissipation. Again, the line was adjusted to account for noncoherent electron scattering. During dissipation the peak flux of the line profile is reduced from ∼eleven to nine times the continuum flux, but does not go into absorption. Dissipating for an additional 1000 yr further decreases the peak flux of the line profile by a negligible ∼0.1%. Figure 12 shows the Hα EW dissipates too slowly to match observations in this scenario (red line).

The goal of this work is to use VDD theory to model the quasi-steady state of 66 Oph's disk, and then follow the subsequent 21 yr dissipation until the Hα line transitions from emission to absorption. Using thousands of models constrained by photometric, polarimetric, and Hα observations, we determined the physical properties of the disk as it evolved.

We determined a set of stellar parameters (Table 2) using MCMC fitting of the UV spectrum. These parameters are consistent with those reported in the literature and were also confirmed by comparison to observations of 66 Oph obtained after dissipation (see Figure 10).

We used diskless polarimetric observations obtained in 2016 to determine the interstellar polarization in the direction of 66 Oph more accurately than previous studies using nearby field stars (see the work by Draper et al. 2014). We modeled the interstellar polarization by fitting a Serkowski law to these diskless observations. The interstellar polarization was then subtracted from HPOL observations, acquired when the disk was present, to determine the polarization level and position angle of the disk (Figure 16). Using the QU diagram (Figure 2) we also determined the disk's polarization position angle to be θ_int = 98° ± 3°. Our two methods of determining the polarization position angle were found to agree. This confirms our interstellar polarization correction using the 2016 LNA/OPD observations.

In this investigation, we began by studying in detail the quasi-steady-state phase of the disk, prior to the dissipation that started in 1989. Then, we used four different hydrodynamic scenarios to explore disk evolution. These scenarios assumed the disk was either isothermal or nonisothermal, and built to the quasi-steady state before dissipating or dissipated from a defined density distribution.

Our best-fit disk model to the quasi-steady state, constrained by observations of the Hα line profile and SED, and verified by the polarization, has a density of ρ₀ = 2.5 × 10⁻¹¹ g cm⁻³, n = 2.6, and i = 575, with the outer disk radius at the lower limit of R_out = 100 R_eq. The density slope, n, is consistent with the results of Waters et al. (1987), who reported n = 2.5, and Vieira et al. (2017), who found n = 2.6; however, both of these studies found disk base densities at the innermost disk approximately one order of magnitude less dense. The close agreement of the density slope suggests it has been well constrained; however, further study into the base density is necessary.

The scenario assuming an isothermal disk that builds and then dissipates (the constant α T, n = 3.5 scenario) reached a density slope of only n = 3.5, much steeper than the best-fit quasi-steady state model with n = 2.6. As a result, the disk was not bright enough and dissipated too quickly due to the rapid density fall-off with radial distance. Interestingly, this scenario predicted a brightening of the Hα flux shortly after the onset of dissipation. We attribute this to the heating of cool regions of the disk at greater radial distance from the star as the innermost disk accretes. This phenomenon was not observed in disks with smaller values of n.

Next, we explored isothermal disks that dissipated from our best-fit quasi-steady state (the constant α T, n = 2.6 scenario). We were able to constrain R_out simultaneously with α as shown in Figure 13. We found R_out = 115 R_eq and α = 0.4. This scenario and these values successfully reproduced the rate of dissipation for all observables considered (Hα, and V-band polarization and magnitude). A remarkable fit of Hα line profile was obtained at selected epochs of disk dissipation (Figure 14).

Krtička et al. (2011) estimated the outer radius of Be star disks. Using their Equation (5) and estimating the rate of change in the moment of inertia using the prescription from Claret & Gimenez (1989) along with our determination of $\dot{M}$ from SINGLEBE, we obtained an outer disk radius of ∼125 R_eq, which is in relatively good agreement with our value determined above.

For the scenario with a nonisothermal disk that builds and dissipates (the variable α T, n = 2.7 scenario) we used a power-law α T that varies with radius as r⁻². The viscous timescales of this scenario were found to be too large to reproduce the rate of dissipation. After 21 years of dissipation, much of the disk remained, and while the visible and near-IR flux matched the diskless star, the Hα and far-IR remained too bright. However, the outer disk dissipated very little, and the predicted radio emission closely matched the JVLA observations from 2010.

Finally, we considered a nonisothermal disk and followed dissipation again from our best-fit quasi-steady state (the variable α T, n = 2.6 scenario). We used the same procedure as the constant α T, n = 2.6 scenario, and the power-law α T from the variable α T, n = 2.7 scenario. This scenario dissipated similar to the other variable α T scenario, and the predicted observables revealed that the disk was still present after 21 years.

Our models showed that the excess radio emission observed in 2010 (Figure 6) can only be reproduced if 66 Oph's outermost disk is unaffected by the dissipation. This suggests that a large mass reservoir was built in the outer disk, as discussed in Rímulo et al. (2018) and Ghoreyshi et al. (2018).

Overall, we found that the constant α T, n = 2.6 (isothermal) scenario, using α = 0.4 and R_out = 115 R_eq, best reproduces the observed dissipation of 66 Oph's disk. However, this scenario fails to explain how the disk could build to n = 2.6. These smaller values of n can be reached with a nonisothermal assumption; however, all of our nonisothermal scenarios failed to match the timescale of dissipation. It is possible that prescriptions for α T other than a power law might help this issue.

Another explanation for low values of n might be the accumulation effect described by Panoglou et al. (2016) and by Cyr et al. (2017). Basically, the disk density becomes shallower due to the build-up of material inward of a binary companion. This effect was not explored here, but we point out that our results indicate that the disk around 66 Oph is quite large, suggesting that the existence of a close by companion is unlikely.

In future work, we aim to better understand the structures of α and T in the disk. It may be possible to empirically model α as a function of radius using a set of other time-series emission lines during the dissipation phase, as different lines form within different radial locations in the disk.

The authors would like to thank the anonymous referee for careful reading and helpful comments that have greatly improved this paper. C.E.J. wishes to acknowledge support though NSERC, the Natural Sciences and Engineering Research Council of Canada. This work was made possible through the use of the Shared Hierarchical Academic Research Computing Network (SHARCNET). A.C.C. acknowledges support from CNPq (grant 311446/2019-1) and FAPESP (grant 2018/04055-8). Development and testing of the routines used was completed on 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. A.C.R. acknowledges the support of FAPESP grants 2017/08001-7 and 2018/13285-7. This work makes use of observations from the LCOGT network. K.C.M. would like to thank Daniel Bednarski and Matheus Genaro from the Instituto de Astronomia, Geofísica e Ciências Atmosféricas, at the Universidade de São Paulo, for reducing the Hα and polarization observations made at OPD, and guidance in working with the emcee code, respectively. The authors would also like to thank Thomas Rivinius from the European Southern Observatory for providing the set of hydrogen spectral line observations used in this work, Chris Tycner from Central Michigan University for a set of Hα equivalent width measurements, and Geraldine Peters for her generosity and willingness to share a set of Hα observations.

Software: emcee (Foreman-Mackey et al. 2013), HDUST(Carciofi & Bjorkman 2006), SINGLEBE (Okazaki 2007).