© E. S. G. de Almeida et al. 2020

1 Introduction

Classical Be stars are main-sequence B-type stars that show (or showed at some time) Balmer lines in emission and infrared excess in their spectral energy distribution. The Be phenomenon is found among the entire spectral range of B stars (e.g., Townsend et al. 2004): M_⋆ from ~3 M_⊙ (B9, T_eff ~ 12 000 K), up to ~18 M_⊙ (B0, T_eff ~ 30 000 K). These observational characteristics are well explained as arising from a dust-free gaseous disk that is supported by rotation with a slow radial velocity (see, e.g., Rivinius et al. 2013). The most successful theory to explain the evolution of the disk structure is the so-called viscous decretion disk (VDD) model, where its dynamics are driven by viscosity (e.g., Lee et al. 1991; Okazaki 2001; Bjorkman & Carciofi 2005).

It is widely accepted that fast rotation plays an important role in the formation of the Be star disk. However, while interferometric analyses typically provide rotational rates v_rot ∕v_crit ≳ 0.7 (e.g., Meilland et al. 2012; Cochetti et al. 2019), some statistical studies show rates ranging from ~0.3 up to 1.0 (e.g., Cranmer 2005; Zorec et al. 2016). Moreover, it is still not clear whether the rotational rate is correlated to other stellar parameters such as the effective temperature (e.g., Cochetti et al. 2019). Hence, despite the success of the VDD model, the physical mechanism(s) driving the mass injection remains unclear and a detailed physical characterization for the central star and the disk structure is mandatory to better understand the Be phenomenon. By gaining access to geometry on the milliarcsecond scale and kinematics on a few tens of km s⁻¹ scale, spectro-interferometry offers a unique opportunity to probe the circumstellar environment and stellar surfaces of Be stars (see, e.g., Chesneau et al. 2012; Stee & Meilland 2012).

The bright, late-type Be star (type B7IVe) o Aquarii (HD 209409) is known to have a fairly stable disk (Sigut et al. 2015). The stability of the circumstellar disk is evidenced by the quasi-constant equivalent width in the Hα line, double-peak separation, and the absence of long-term violet-to-red (V∕R) peak variations (e.g., Rivinius et al. 2006; Sigut et al. 2015). This star shows a high value of v sin i ~ 282 km s⁻¹ (Frémat et al. 2005) and a shell absorption in Hα, thus indicating a high stellar inclination angle of about 70°, as discussed below.

Meilland et al. (2012) presented the first spectro-interferometric analysis of o Aquarii with the VLTI/AMBER instrument as part of their AMBER survey of eight bright Be stars. Despite the low data quality and very limited number of observations (just one measurement), they were able to significantly constrain the disk geometry and kinematics. They found that the disk emission in the Brγ line, modeled as an elliptical Gaussian distribution,had a FWHM of 14 ± 1 D_⋆ (with R_⋆ = 4.4 R_⊙), where D_⋆ and R_⋆ are, respectively, the stellar diameter and radius. They estimated the inclination angle as i = 70 ± 20° and found a stellar rotational rate of v_rot ∕v_crit = 0.77 ± 0.21 (Ω/Ω_c = 0.93${}_{-0.17}^{+0.06}$), where v_crit and Ω_crit are, respectively,the linear and angular critical velocity. New VLTI/AMBER spectro-interferometric measurements of o Aquarii were presented in the Be star survey of Cochetti et al. (2019). Here, they obtained seven good-quality measurements for o Aquarii (i.e., 21 baselines). Using a similar model as in Meilland et al. (2012), they found a Brγ emission FWHM significantly smaller than in Meilland et al. (2012), 8 ± 0.5 D_⋆ (with R_⋆ = 4.4 R_⊙), and better constrained the object inclination angle (70 ± 5°).

A detailed analysis of o Aquarii using Hα spectroscopy and interferometry was performed by Sigut et al. (2015). These authors combined large band (15 nm) interferometric data centered on Hα, obtained from the Navy Precision Optical Interferometer (NPOI), with Hα spectroscopy from the Lowell Observatory Solar-Stellar Spectrograph. Using the radiative transfer code BEDISK (Sigut & Jones 2007), they were able to reproduce simultaneously the visibility, Hα line profile, and spectral energy distribution (SED), and showed that the disk is quite stable for up to about ten years. Interestingly, they found a disk extension in Hα (Gaussian FWHM of 12.0 ± 0.5 D_⋆) close to the one determined by Meilland et al. (2012) in Brγ (FWHM of 14 ± 1 D_⋆). They concluded that this is uncommon since most previously studied Be stars exhibit a larger (up to two times) disk emission region in Hα than in Brγ.

In this paper, we present new CHARA/VEGA spectro-interferometric measurements of o Aquarii centered on the Hα emission line (λ = 0.656 μm). They are analyzed conjointly with the AMBER Brγ line (λ = 2.166 μm) measurements from Meilland et al. (2012) and Cochetti et al. (2019), using models of increasing complexity: simple geometric models, kinematic models, and radiative transfer models. This is the first time the code HDUST has been used to model simultaneously spectro-interferometric data from Hα and Brγ. It is the second time for the kinematic model (i.e., after the δ Scorpii data published in Meilland et al. 2011). This multi-wavelength and multi-line approach allows us to draw a more complete picture of the stellar surface and circumstellar environment of the Be star o Aquarii.

This paper is organized as follows. In Sect. 2, we present the observations and the data reduction process. Our analysis using geometric models of the VEGA calibrated (absolute) visibility is shown in Sect. 3. In Sect. 4, we fit the VEGA and AMBER differential visibility and phase with a kinematic model using a Markov chain Monte Carlo (MCMC) model fitting method. In Sect. 5, all the interferometric data are analyzed in terms of 3D non-LTE radiative transfer models. Our kinematic and radiative transfer models are discussed in Sect. 6. In Sect. 7, our best-fit models are compared to non-interferometric observables: the spectral energy distribution and line profiles (Hα and Brγ). The comparison with polarimetric data is performed in Sect. 8.4.3 in the context of the disk stability. In Sect. 8, we discuss the morphological, kinematic, and physical descriptions for o Aquarii and its circumstellar disk. Our conclusions are summarized in Sect. 9.

2 Observations

2.1 CHARA/VEGA

The VEGA instrument (Mourard et al. 2009) is one of the two visible beam combiners on the CHARA Array (ten Brummelaar et al. 2005). It can simultaneously combine up to four beams, operating at different wavelengths from 450 to 850 nm. VEGA is equipped with two cameras (blue and red detectors) that can observe in two different spectral domains simultaneously (around the Hβ and the Hα lines). Currently, it is the only instrument at the CHARA Array with a spectral resolution high enough to resolve narrow spectral features such as atomic and molecular lines. It offers 3 spectral modes: R = 1000 (LR), R = 6000 (MR), and R = 30 000 (HR).

o Aquarii was observed 50 times with VEGA between 2012 and 2016 in MR mode centered on the Hα emission line at 0.656 μm. The 2012 and 2016 observations were focused on the disk geometry and kinematics and data were taken with small baselines (up to 105 m) and without stellar calibrators. On the other hand, the 2013 and 2014 campaigns were aimed at constraining not only the Hα emission, but also the R-band continuum geometry. Consequently, observations were carried out with longer baselines (up to 330 m) with a standard calibration plan alternating observations of the science target and few calibrator stars chosen using the SearchCal (Bonneau et al. 2006) tool developed by the Jean-Marie Mariotti Center (JMMC)¹. Table A.1 shows useful information about the stars used as interferometric calibrators during these campaigns. The complete log of observations is presented in Table A.2 and the corresponding uv plane coverage for the VEGA observations is plotted in Fig. 1.

Data were reduced using the standard VEGA data reduction software² described in Mourard et al. (2012). For all programs, differential visibility and phases were computed from the intercorrelation between a fixed 15 nm window centered on Hα and a sliding smaller window (i.e. 1, 2, or 5 Å, depending on the data quality). For the 2013 and 2014 data, the raw squared visibility was computed for o Aquarii, and its calibrators, using the auto-correlation method on a 15 nm band centered on the Hα emission line (649–664 nm) and another band in the close-by continuum (635–650 nm). Then the transfer function was estimated assuming the diameter of the calibrators recorded before and after the science target observation, and its uncertainty using a weighted standard deviation. Finally, for each measurement, the calibrated squared visibility was derived by dividing o Aquarii’s raw squaredvisibility by the estimated transfer function.

Fig. 1 uv plan coverage obtained around Hα (0.656 μm) with CHARA/VEGA (green) and Brγ (2.166 μm) with VLTI/AMBER (red).

2.2 VLTI/AMBER

The AMBER instrument (Petrov et al. 2007) was a three-beam combiner (decommissioned in 2018) at the Very Large Telescope Interferometer (VLTI). It operated in the H- and K-bands with three spectral resolutions: R = 35 (LR), R = 1500 (MR), and R = 12 000 (HR). It offered the highest spectral resolution at the VLTI, being the most adapted for studying the gaseous environment in emission lines.

o Aquarii was observed with AMBER during two observing surveys of Be stars in 2011 (ESO program 087.D-0311) and in 2014 (ESO program 094.D-0140). The observations were performed in HR mode in K-band centered on the Brγ emission line at 2.166 μm. The data from 2011 was published in Meilland et al. (2012) and the 2014 data in Cochetti et al. (2019). During this second survey, seven measurements were acquired for o Aquarii with three different triplets. The log of AMBER observations is also presented in Table A.2 and the corresponding uv plane coverage is plotted in Fig. 1.

Calibration was performed using similar methods as the one described for VEGA. However, AMBER measurements were often affected by a highly variable transfer function mainly due to the variable quality of the fringe tracking performed by the FINITO fringe tracker, during the long exposure time needed to perform HR mode observations. As it was the case during our o Aquarii observations, we present in this paper only the analysis of differential measurements obtained using the standard AMBER data reduction software amdlib (Tatulli et al. 2007; Chelli et al. 2009).

3 Geometric modeling: VEGA calibrated visibility

In this section, we fit the Hα and continuum squared visibilities (V²) from the VEGA observations where calibrators were observed. We note that as the AMBER data were not calibrated, such analysis cannot be performed on the K-band continuum and Brγ line.

To determine if we can separate the circumstellar disk and the stellar photosphere emissions and constrain their geometry independently, we fitted our data with geometric models of increasing complexity: one-component models (uniform disk, UD, or a uniform ellipse) and two-component models (UD plus UD, Gaussian disk, or uniform or Gaussian ellipse).

Here, the first component represents the stellar surface and the second one the circumstellar disk. To perform our fit, we used the LITpro model fitting software (Tallon-Bosc et al. 2008) for optical and infrared interferometric observations developed by the Jean-Marie Mariotti Center (JMMC)³.

In Fig. 2, we show the comparison between the visibility curves of our best-fit models to the VEGA data both in the continuum and Hα bands. One sees that the object is partially resolved in the continuum and the Hα line. The lower level of the visibility in the band centered on the Hα line clearly shows that the object is larger in Hα than in the close-by continuum region. Assuming that the emission originates from both the stellar photosphere and a circumstellar disk, the lower visibility in Hα is due to a larger fraction of the Hα flux coming from the disk than from the star. In contrast, the flux contribution from the star is greater than that from the disk in the continuumR-band.

Our main results are summarized in Table 1. We only show our results using UD models since there is no improvement in terms of reduced χ² (${\chi }_{\text{r}}^2$) when considering more complex models, that is, with a higher number of free parameters. For the continuum band, there is no significant improvement in terms of reduced χ² between a simple UD and a two-component UD model. The central star is clearly resolved by the longer baselines and its extension is significantly constrained with a UD diameter of θ = 0.28 ± 0.01 mas (${\chi }_{\text{r}}^2$ ~ 1.1). This value corresponds to an upper limit to the stellar diameter measurements neglecting the putative contribution of the circumstellar disk in the R-band continuum. Adding a second component to the model only marginally reduces the extension of the first component. The contribution of the second component, representing the circumstellar disk, is small (F₂ = 0.03 ± 0.03), thus the extension of the disk cannot be constrained.

Unlike the continuum case, the situation is quite different in the band centered on the Hα line. The single uniform disk gives a significantly higher ${\chi }_{\text{r}}^2$ ~ 2.8 for a best-fit model with θ = 0.36 mas. In this case, adding a second component reduces ${\chi }_{\text{r}}^2$ by a factor of two, leading to ${\chi }_{\text{r}}^2~1.3$. Using a model with two uniform disks, we converge to a diameter of the first component similar to the one found from the continuum, that is, 0.26 ± 0.02 mas. The flux contribution of the second component and its extension are significantly constrained. However, the uncertainty remains quite large, that is, F₂ = 0.15 ± 0.03 and θ₂ = 6.5 ± 2.1 mas (see Table 1).

Considering that the first component of our model represents the stellar photosphere, our measurement is slightly higher than the value assumed in the work of Sigut et al. (2015) of 0.22 mas. However, their adoption for the stellar angular diameter is based on a spectral type-radius relation for B dwarf stars (Townsend et al. 2004). Moreover, this value of 0.22 mas represents the polar radius. o Aquarii is a fast rotator likely to be significantly flattened, and our measurements are spread over different orientations, so that we end up measuring a mean radius of the star projected on the sky. Assuming a distance of 144 pc (derived from the Gaia DR2 parallaxes, Gaia Collaboration 2018), θ₁ = 0.26 ± 0.02 mas corresponds to a stellar radius R_⋆ = 4.0 ± 0.3 R_⊙.

Finally, to try to detect any possible stellar or circumstellar disk flattening from the squared visibility measurements, we also computed individual uniform disk equivalent diameter for each V² measurement. This analysis of the uniform disk diameter for o Aquarii, as a function of the VEGA baseline orientation, is shown in Fig. 3. As expected from our analysis (considering uniform elliptical models), we do not find any evidences of flattening from modeling our V² dataset since no clear trends are found in the model residual as varying the baseline position angle.

Fig. 2 VEGA V2 measurements in the close-by continuum band (top) and in the Hα band (bottom) are shown in red points. Our best-fit models consisting of one (solid line) and two (dashed line) uniform disks are overplotted in blue. See Table 1 and text for discussion.

Fig. 3 Top panels: uniform disk diameter derived from each individual VEGA V2 measurements (continuum band in the left and Hα band in the right) plotted as a function of the baseline position angle (PA). The red dotted line represents the best-fit diameter from modeling all the data in each band (θ = 0.28 mas in the continuum and θ = 0.36 mas in the Hα band). Bottom panels: corresponding normalized residuals.

4 Kinematic modeling: VEGA and AMBER differential data

To constrain the geometry and kinematics of the circumstellar gas in the Hα and Brγ lines, we fit the VEGA and AMBER differential visibility and phase measurements using a simple bi-dimensional kinematic model for a rotating disk⁴.

4.1 The kinematic model

This kinematic model was already used in a series of papers about spectro-interferometric modeling of Be stars, including Delaa et al. (2011), Meilland et al. (2012), and Cochetti et al. (2019), and is presented in detail in these references.

In short, the intensity map for the central star is modeled as a uniform disk, and the circumstellar disk as two elliptical Gaussian distributions, one for the flux in continuum, and the other one for the flux in line. The disk is geometrically thin so that the ellipse flattening ratio is set to 1∕cosi, where i is the inclination angle. The disk intensity map in the line is computed taking into account the Doppler effect due to the disk rotational velocity in the considered spectral channels. The parameters of our kinematic model are the following:

• (i)

  The simulation parameters: size in pixels (n_xy), field of view in stellar diameters (fov), number of wavelength points (n_λ), central wavelength of the emission line (λ₀), step size in wavelength (δλ), and spectral resolution (Δλ).

• (ii)

  The global geometric parameters: stellar radius (R_⋆), distance (d), inclination angle (i), and disk major-axis position angle (PA).

• (iii)

  The disk continuum parameters: disk major-axis FWHM in the continuum (a_c), disk continuum flux normalized by the total continuum flux (F_c).

• (iv)

  The disk emission line parameters: disk major-axis FWHM in the line (a_line) and line equivalent width (EW).

• (v)

  The kinematic parameters: rotational velocity (v_rot) at 1.5 R_p (polar radius) and exponent of the rotational velocity power-law (β).

4.2 Model fitting using the MCMC method

To perform our model fitting, we used the code emcee (Foreman-Mackey et al. 2013). This is an implementation in Python of the MCMC method from Goodman & Weare (2010). Some recent works on stellar interferometry used this code (see., e.g., Monnier et al. 2012; Domiciano de Souza et al. 2014, 2018; Sanchez-Bermudez et al. 2017).

The simulation parameters were set as follows: n_xy = 256, fov = 60 D_⋆, n_λ = 60 (VEGA) and 110 (AMBER), λ₀ = 6563 Å (VEGA) and 21 661 Å(AMBER), δλ = 2.5 Å (VEGA) and 1.0 Å (AMBER), and Δλ = 5.0 Å(VEGA) and 1.8 Å (AMBER). To reduce the number of free parameters, we set R_⋆ = 4.0 R_⊙ and d = 144 pc. We also fixed the disk continuum extension a_c and flux F_c to 0 for VEGA (i.e., neglecting the disk contribution in the continuum, based on our analysis of the VEGA V² data). In the AMBER analysis, we adopted a_c = 3 D_⋆ and F_c = 0.2 from Cochetti et al. (2019). The line equivalent width was set to 19.9 Å in Hα (Sigut et al. 2015). For Brγ, we computed the EW using the AMBER spectra from all observations and found a mean value of 13.6 ± 1.1 Å, which is compatible with the value from Meilland et al. (2012), 12.6 Å, but not with the result from Cochetti et al. (2019) of 18.1 Å. Finally, from the ten parameters of the kinematic model, the fitting of the VEGA and AMBER data were performed with at most five free parameters: i, PA, a_line, v_rot, and β.

The likelihood function (p_like) of the MCMC procedure was chosen as $\ln (p_{\text{like}})=-{\chi }_{\text{total}}^2/2$, where ${\chi }_{\text{total}}^2$ is the sum of the χ² computed for the differential visibility and the differential phase. Thus, our attempt to converge to samples of parameters that maximizes the likelihood function means the minimization of the total χ² between our interferometric data and the kinematic model.

We performed three different model fitting tests with different constraints on the value of v_rot:

• (i)

  Five free parameters: i, PA, a_line, v_rot, and β. Without theinclusion of any prior probability function in the analysis.

• (ii)

  Fourfree parameters: i, PA, a_line, and β. The stellar rotational velocity v_rot is fixed on the critical value of 391 km s⁻¹ (Frémat et al. 2005).

• (iii)

  Five free parameters: i, PA, a_line, v_rot, and β. We take into account a prior probability function p_prior on v sin i. Adopting μ = 282 km s⁻¹ and σ = 20 km s⁻¹, from the measured v sin i = 282 ± 20 km s⁻¹ (Frémat et al. 2005), we have the following expression for p_prior: $$\text{ln}(p_{\text{prior}})=\frac{-{(v\sin i-\mu )}^2}{2{\sigma }^2},$$(1)

  where v sin i is calculated from the sampled MCMC values for the stellar rotational velocity and inclination angle.

Hence, considering a high weight on p_prior, the following quantity for the posterior probability function p_post is maxi- mized: $$\ln (p_{\text{post}})=-100\left(\frac{{(v\sin i-\mu )}^2}{2{\sigma }^2}\right)-\frac{{\chi }^2}{2}.$$(2)

Note that this is equivalent to the case of equal weights for p_prior and p_like, but considering a lower error bar on v sin i, namely, σ = 2 km s⁻¹.

We typically used several hundreds of walkers (~300–900) for the MCMC run. Convergence was obtained for about 50 to 100 iteration steps in each walker, but we used a conservative value of 150 steps in the burn-in phase and 50 in the main phase to estimate the parameters values and uncertainties. Overall, we found a mean acceptance fraction of ~0.5–0.6 in our MCMC tests. This is close to the optimal range for this parameter of ~0.2–0.5 (see, e.g., Foreman-Mackey et al. 2013).

4.3 Best-fits in Hα and Brγ

We modeleda total of 117 (VEGA) and 24 (AMBER) measurements of differential visibility and phase. The best-fit parameters for the MCMC fit with a prior on v sin i (test iii, described above) are presented in Table 2. The corresponding histograms and the two-by-two parameter correlations from this MCMC run (one for VEGA and other for AMBER) are shown in Fig. 4. The corresponding histograms and correlation plots for the other two fits (tests i and ii) are shown in Figs. B.1 and B.2. One sees that the values of i, PA, and a_line, derived from each emission line, differ only marginally in all the fitting tests, showing the robustness of the solution for these parameters.

In Fig. 5, we show examples of VEGA and AMBER data in comparison to our best-fit kinematic models. For later discussion in Sect. 6, the visibility and phase from our best-fit HDUST model is also presented here. Our best-fit kinematic models are able to reproduce both the VEGA and AMBER differential data well. We found a reduced χ² of ~4.0 and 1.6 from fitting, in a separate way, respectively, the VEGA and AMBER datasets.

We derived compatible values for the disk PA (~110°) from fittingthe VEGA and AMBER data with an uncertainty up to ~2°. This result agrees well with previous studies (e.g., Meilland et al. 2012; Touhami et al. 2013; Sigut et al. 2015; Cochetti et al. 2019). On the other hand, the inclination angle determined from the fit to the VEGA data is significantly smaller (i = 61.2 ± 1.8°) in comparison to the one determined from fitting AMBER (i = 75.9 ± 0.4°). This latter value is in good agreement with the results for i found by Meilland et al. (2012) and Cochetti et al. (2019). We also constrain the disk extension with a good precision: a_line = 10.5 ± 0.3 D_⋆ in the Hα line and a_line = 11.5 ± 0.1 D_⋆ in the Brγ line. These values are compatible with the ones determined by Sigut et al. (2015) in Hα and Meilland et al. (2012) in Brγ.

Another aspect concerning the disk extension in Brγ is the significant discrepancy seen in comparison to a_line = 8.0 ± 0.5 D_⋆ from Cochetti et al. (2019). However, these authors used a larger value for the stellar radius of 4.4 R_⊙ and a closer distance of 134 pc (van Leeuwen 2007), having thus the angular size of the stellar diameter larger in ~19% than the one assumed in our kinematic analysis from our results in Sect. 3. Considering all the other parameters fixed, this results in a smaller disk extension in ~19% than one found from our analysis. Nevertheless, the largest contribution to this discrepancy between our results and the ones from Cochetti et al. (2019) is due to their high value of equivalent width in the Brγ line of 18.1 Å, as discussed in Sect. 4.2, that also implies in a smaller disk extension in this line.

From our various tests, we showed that β and v_rot are strongly correlated. To precisely determine their dependence, we computed a grid of kinematic models varying just these two parameters in a regular step size. The values for i, PA, a_line are fixed from Table 2. The resulting ${\chi }_{\text{r}}^2$ maps are shown in Fig. 6. As expected, one sees that v_rot and β are highly correlated for the VEGA and AMBER data. This high degeneracy can be understood since these two parameters provide the rotational velocity structure in the disk: it is hard to distinguish the effects of each one on the modeling of spectro-interferometric (and spectroscopic) data.

Furthermore, we see that β = −0.5 (Keplerian disk) provides unrealistically high values for the stellar rotational velocity ((≳)400 km s⁻¹ ; gray region) of o Aquarii (VEGA analysis). For AMBER, v_rot is significantly reduced to about 300–400 km s⁻¹. As shown in Fig. 6, our results from AMBER are consistent with a nearly Keplerian rotating disk (β ~ 0.43). However, it is conspicuous that the β value calculated from the VEGA data (β ~ 0.30) shows such a large departure from the Keplerian case.

Cochetti et al. (2019) derived a stellar rotational velocity of 355 ± 50 km s⁻¹ and β = −0.45 ± 0.03. This is in fair agreement with our results for both v_rot and β. Considering our MCMC test (ii), where v_rot is fixed to the critical value and β is a free parameter, the results for β are shifted to higher values (more positive) with β ~ −0.42 (VEGA) and −0.54 (AMBER).

Therefore, regardless the MCMC fitting considered here, we verify a discrepancy of about 0.1 between the value of β derived from the Hα and Brγ lines. Our results from the AMBER analysis (Brγ) seems to be consistent with a nearly Keplerian rotating disk, but we verified a larger departure from β = −0.5 for the VEGA analysis (Hα).

Fig. 4 Histogram distributions and two-by-two correlations (after the burn-in phase) for the free parameters of our best-fit kinematic models using MCMC for the VEGA (left panel) and AMBER (right panel) differential data. The median values are shown in solid red lines and the first and third quartiles in dashed red lines. The median and the first and third quartiles estimated for the parameters of our best-fit models (VEGA and AMBER) are presented inTable 2. In the correlation plots, darker points correspond to models with lower values of χ2. See text for discussion.

Fig. 5 Comparison between our best-fit kinematic models (dashed red; Table 2) and two different VEGA (top panels) and AMBER (bottom panels) measurements (black line). Our best-fit HDUST model is also shown (dashed blue; Table 5; discussion in Sect. 6). δλ of the kinematic model and AMBER data is increased to 1.8 Å in order to compare them to the HDUST model (δλ fixed to 1.8 Å).

Fig. 6 ${\chi }_{\text{r}}^2$ maps of 40 000 kinematic models as a function of vrot and β from the fit to VEGA (top panel) and AMBER (bottom panel) differential data. Only these two parameters were varied in a regular step in the intervals shown here. The other parameters are fixed (Table 2). Our results found from the MCMC analysis for vrot and β are indicated with red crosses. In order to highlight the correlation between β and vrot, the gray region corresponds to an arbitrary number of models, encompassing about the 5000 best models in both cases. The value of β = −0.5 (Keplerian disk) and our determination for vrot are marked in dashed black line. Note the strong correlation between the stellar rotational velocity and the disk velocity law exponent in both the cases. Also, note that a Keplerian disk is found from modeling the AMBER data, but not from VEGA.

5 Radiative transfer modeling

5.1 The code HDUST

We used the 3D non-LTE radiative transfer code HDUST⁵ (Carciofi & Bjorkman 2006, 2008) to perform a deeper physical analysis of o Aquarii. In addition to geometric and kinematic parameters, we seek to derive the density and temperature distributions in the disk, and the spectral energy distribution (SED), none of which was provided by the two simpler models considered in the two previous sections. HDUST uses a Monte Carlo method to solve the radiative transfer, statistical and radiative equilibrium equations for arbitrary density and velocity distributions in gaseous (pure hydrogen) or dusty circumstellar environments.

This code is well-suited to model the circumstellar environment of Be stars as it implements the VDD model. Thus, the disk velocity law is assumed to be Keplerian (β fixed to −0.5). Many previous studies explored formal solutions of the VDD model in several limiting cases. For example, Bjorkman & Carciofi (2005) investigated the isothermal, steady-state case of a disk formed by a steady mass injection rate over a long time. Effects due to non-isothermal temperature structure were studied by Carciofi & Bjorkman (2008). Haubois et al. (2012) studied the temporal evolution of the disk structure that is subject to variable mass inject rates. Finally, the effects of a binary companion on the disk were studied by Okazaki et al. (2002), Oudmaijer & Parr (2010), Panoglou et al. (2016), and Cyr et al. (2017), among others.

From these studies, the radial density profile in Be star disks is found to be quite complex, for example, depending on the disk age, dynamical state, or presence of a binary companion. Despite this complexity, several studies have shown that the global behavior of this density profile is successfully approximated by a simple radial power-law (e.g., Touhami et al. 2009; Vieira et al. 2017). Considering also that the vertical density structure is that of an isothermal disk (hydrostatic assumption in the z-axis), the disk density can be parameterized as follows: $$\rho (r,z)={\rho }_0{\left(\frac{R_{\text{eq}}}{r}\right)}^m\exp \left(\frac{-z^2}{2H{(r)}^2}\right),$$(3)

where ρ₀ is the disk base density, R_eq is the equatorial radius, and H(r) is the (isothermal) disk scale height given by: $$H(r)=H_0{\left(\frac{r}{R_{\text{eq}}}\right)}^{3/2},$$(4)

and H₀ is the scale height at the disk base, $$H_0=c_{\text{s}}{R}_{\text{eq}}{\left(\frac{G{M}_{\star }}{R_{\text{eq}}}\right)}^{-1/2},$$(5)

where M_⋆ is the stellar mass, G the gravitational constant, and c_s the sound speed velocity which depends on the local disk temperature T: $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu m_{\text{H}}}},$$(6)

where k_B is the Boltzmann constant, μ is the mean molecular weight of the gas, m_H is the hydrogen mass, and T is adopted as 0.72 T_pol, where T_pol is the polar effective temperature (see Correia Mota 2019).

HDUST has been used a few times to model spectro-interferometric observations (e.g., Carciofi et al. 2009; Klement et al. 2015; Faes 2015). From the solution of the radiative transfer problem, we are able to calculate synthetic spectra and intensity maps as a function of the wavelength around specific spectral lines. We estimated the stellar and circumstellar disk parameters from the comparison of our spectro-interferometric observations (visible and near-infrared) with synthetic observables computed from the Fourier transform of HDUST monochromatic intensity maps.

5.2 BeAtlas grid

Since a few hours are needed to compute a single HDUST model, it is not possible to perform an iterative model fitting procedure similar to the one described in Sect. 4. To overcome this issue, we used a pre-computed grid of HDUST models called BeAtlas (Faes 2015; Correia Mota 2019). The BeAtlas grid is presented and described in detail by these references. It consists of ~14 000 models with images (specific intensity maps), SEDs, and spectra calculated in natural and polarized spectra, over several spectral regions, including the Hα and Brγ lines that are of interest for the analysis of our VEGA and AMBER dataset.

In Table 3, we show the parameter space covered by BeAtlas. Five physical parameters are varying in the grid. The stellar mass M_⋆, the inclination angle i, and the stellar oblateness R_eq∕R_p, fully describe the star. Other stellar parameters such as the stellar polar radius (R_p), rotational velocity (v_rot) and linear and angular rotational rates (v_rot∕v_crit and Ω∕Ω_crit) can be computed from M_⋆ and R_eq∕R_p assuming rigid rotation under the Roche model (see, e.g., Carciofi & Bjorkman 2008). The two last parameters in Table 3describe the circumstellar disk structure and are parameterizations of the VDD model: the base surface density (Σ₀) and the radial density exponent (m).

The previously described volume mass density (Eq. (3)) and the surface mass density are related as follows: $$\Sigma (r)\equiv {\displaystyle {\int }_{-\infty }^{+\infty }\rho }(r,z)\text{d}z,$$(7) $$\rho (r,z)=\frac{\Sigma (r)}{H(r)\sqrt{2\pi }}\exp \left(\frac{-z^2}{2H{(r)}^2}\right).$$(8)

From that, to facilitate the comparison to other disk models, we note that the relation between the volume and surface mass densities at the base of the disk is given by: $${\rho }_0={\Sigma }_0\sqrt{\frac{G{M}_{\star }}{2\pi c_{\text{s}}{}^2{R}_{\text{eq}}{}^3}}.$$(9)

The range of values for Σ₀ and m in the grid encompasses somewhat extreme cases in the literature for the circumstellar disk of Be stars. For example, see Fig. 7 of Vieira et al. (2017). The listed values of Σ₀ correspond to ρ₀ from ~10⁻¹² g cm⁻³ to ~10⁻¹⁰ g cm⁻³. Parametric models with m = 3.5 are equivalent to the steady-state solution of the viscous diffusion equation considering an isothermal disk scale height. Thus, concerning the mass density law exponent m, models with m > 3.5 would represent a disk in an accretion phase, while the ones with m < 3.5 a disk in an ongoing process of dissipation (see, e.g., Haubois et al. 2012; Vieira et al. 2017).

5.3 Results

We performed four different analyses of our data using different subsets. For that, the reduced χ² between the predicted interferometric observables from each HDUST model and the data was calculated as follows:

• (i)

  calibrated VEGA V ² in the 642.5 nm band (close-by continuum to Hα).

• (ii)

  VEGA differential visibility and phase (Hα line).

• (iii)

  AMBER differential visibility and phase (Brγ line).

• (iv)

  All the quantities above analyzed together.

Analysis (i) was performed to evaluate the constraint on the stellar mass M_⋆ and oblateness R_eq ∕R_p. In Fig. 7, we show the lowest value of ${\chi }_{\text{r}}^2$ for each value of stellar oblateness and mass from fitting the VEGA V² data in the continuum band. The predicted V² from our best-fit BeAtlas model (with M_⋆ = 4.2 M_⊙ ; Table 5) is overplotted to the VEGA measurements. For comparison, the predicted visibility curve from the BeAtlas model with the highest stellar mass, M_⋆ = 14.6 M_⊙, is also overplotted to the data. These two models have the same values of i, R_eq ∕R_p, Σ₀, and m. In Sect. 3, we presented a similar analysis, but in terms of simple geometric models. For better visualisation, we show in Fig. 7 the local regression fits of ${\chi }_{\text{r}}^2$ as a function of R_eq∕R_p and M_⋆. Like all such calculations in this paper, all these regression fits of ${\chi }_{\text{r}}^2$ are performed with the LOESS method⁶.

As in the analysis with geometric models, we cannot constrain the stellar oblateness using VEGA V² data. On the other hand, the mass is better constrained with M_⋆ ~ 4.8 M_⊙ (B6 dwarf). From Fig. 7, one sees how the measured V² are mismatched by the HDUST model with M_⋆ = 14.6 M_⊙ (unrealistic mass value for o Aquarii) due to the larger polar radius of ~7.4 R_⊙ in this model. Among all the values for M_⋆ in the grid, M_⋆ = 4.2 M_⊙ corresponds to a B7 dwarf star (Townsend et al. 2004). Since o Aquarii shows luminosity class III-IV, it could be expected to have a mass somewhat higher than a dwarf of same spectral type, which is compatible with our results.

In Fig. 8, we show the lowest ${\chi }_{\text{r}}^2$ for each value of disk major-axis position angle PA from the fit to the VEGA and AMBER differential visibilities and phases: analyses (ii) and (iii). Here, the stellar mass is fixed to M_⋆ = 4.2 M_⊙ from analysis (i), which also allows a better comparison to other studies of o Aquarii (e.g. Sigut et al. 2015). In both cases, ${\chi }_{\text{r}}^2$ of the models is minimized around PA = 110°, a value that we adopt in the remaining of this section. This is in good agreement to our results found with the kinematic model in Sect. 4.

In Fig. 9, we present our results from modeling the VEGA and AMBER differential visibility and phase in a separate way – analyses (ii) and (iii) – as well as from the simultaneous fit to all the interferometric data (analysis (iv)). The lowest ${\chi }_{\text{r}}^2$ is shown as a function of the following HDUST parameters: the inclination angle, stellar oblateness, base disk surface density, and the radial disk density law exponent. In Table 4, we show the statistics from these parameters calculated from the HDUST models within a certain threshold of ${\chi }_{\text{r}}^2$, which, in each case, is chosen to match a similar number of models (~15–20 best-models). In Table 4, the parameters for the models with the lowest value of ${\chi }_{\text{r}}^2$ are also shown. In Table 5, we show the parameters for the best BeAtlas model to explain simultaneously all our different interferometric datasets.

Since our HDUST analysis is limited to the pre-computed BeAtlas grid (limited parameter space and selected parameter values), we stress that the results presented here do not correspond to the real χ² minimum toexplain our datasets in the framework of HDUST. Furthermore, the values for the standard deviation are shown in parenthesis in Table 4 since these are not determinations for the error bars on the parameters. They are just an evaluation for the dispersion on the parameters values of the BeAtlas best-models (within in a certain threshold of ${\chi }_{\text{r}}^2$). For example, from fitting AMBER, we found that all the BeAtlas models have Σ₀ = 0.12 g cm⁻², and m = 3.0, up to, respectively, the top 207% and top 240% best-models. For this reason, it is shown, in this case, null standard deviation in Table 4 for these parameters (top 54% best-models).

From the separate analysis of the VEGA and AMBER differential datasets, we are able to describe the stellar and disk parameters, in Hα and Brγ, by the same HDUST model with: R_eq∕R_p = 1.45, Σ₀ = 0.12 g cm⁻², and m = 3.0. One clear exception is found for the inclination angle. From the Hα analysis, ${\chi }_{\text{r}}^2$ is minimized for i = 56.3°. On the other hand, this is achieved with i = 77.2° in the Brγ line. Such discrepancy of ~20° is in agreement with the one found from our kinematic modeling. As expected, the joint analysis to all the data provides an intermediate mean value of ~65° for the inclination angle, showing a larger dispersion (higher standard deviation) in comparison to the results found from the separate analysis for VEGA and AMBER. One sees that the mean value for stellar oblateness is somewhat decreased, when considering all the datasets. However, in this case, the dispersion is significantly increased (±0.09) when compared to the separate VEGA and AMBER differential fits (±0.05–0.07). This happens due to the inclusion of the calibrated VEGA data in the joint analysis that do not allow us to properly infer this parameter (see, again, Fig. 7).

Fig. 7 Analysis (i): lowest value of reduced χ2 (${\chi }_{\text{r}}^2$) for each value of Req∕Rp (left panel) and stellar mass (middle panel) from the HDUST fit to the VEGA V2 data (642.5 nm band). Local regression fits to ${\chi }_{\text{r}}^2$, as a function of the parameter values, are shown in red line. In the right panel, the predicted visibility from our best-fit HDUST model (red points; Table 5) is compared to the VEGA V2 measurements in the continuum band (black points). The predicted visibility from the HDUST model with the highest mass in the BeAtlas grid (14.6 M⊙, highest ${\chi }_{\text{r}}^2$ in the middle panel) is shown in blue points.

Fig. 8 Lowest value of ${\chi }_{\text{r}}^2$ for each value of disk major-axis position angle from the HDUST fit to the VEGA (top panel, analysis (ii)) and AMBER (bottom panel, analysis (iii)) differential visibility and phase. Local regression fits of ${\chi }_{\text{r}}^2$ as a function of the disk PA are shown as a red line.

Fig. 9 Lowest value of ${\chi }_{\text{r}}^2$ for each value of stellar inclination angle, oblateness, base disk surface density, and disk density law exponent from the HDUST fit to the: VEGA differential data (top, analysis (ii)), AMBER differential data (middle, analysis (iii)), and all the interferometric data considered in this section (bottom panel, analysis (iv)). The stellar mass is fixed to 4.2 M⊙ and disk PA to 110°. Local regression fits to ${\chi }_{\text{r}}^2$, as a function of the parameter values, are shown as a red line. The mean parameter values for the sets of best models (Table 4) are marked in dashed black line. Our best-fit BeAtlas model to fit all the interferometric data is shown in Table 5. See text for discussion.

Fig. 10 Intensity maps of our best-fit HDUST and kinematic models at different wavelengths around the Hα line (first two rows) and the Brγ line (last two rows). Flux/pixel is in arbitrary units with the same scale in Hα and Brγ. The image integrated in wavelength around each of these lines (Δλ = 2.7 nm around Hα and 3.9 nm around Brγ) are shown in the second column.

6 Comparison between kinematic and HDUST best-fit models

In Fig. 5, we compare the synthetic differential visibility and phase from our best-fit kinematic and HDUST models to the actual VEGA and AMBER data for a few baselines. Comparisons to non-interferometric observables (spectral energy distribution and line profiles) are presented in Sect. 7. Our best-fit models are compared to all the AMBER data in Fig. C.1. One sees that our best-fit kinematic models do a better job of reproducing both the VEGA and AMBER data. From the separate kinematic modeling of the VEGA and AMBER differential data, the ${\chi }_{\text{r}}^2$ of the model is lower than with HDUST (BeAtlas grid). Fixing the stellar mass to a reliable value for o Aquarii (4.2 M_⊙), our best-fit HDUST model has ${\chi }_{\text{r}}^2$ ~ 6.1 and 4.7 for VEGA and AMBER, respectively. From the kinematic modeling, we found ${\chi }_{\text{r}}^2$ ~ 4.0 and 1.6 to explain these same datasets.

For VEGA, in particular, our best-fit HDUST model adjustment for the measured visibility width is worse than with the kinematic model. This particular issue in modeling the VEGA data can be explained; in HDUST, the disk velocity law exponent is fixed by β = −0.5 (Keplerian disk rotation), while in the kinematic model it is a free parameter. As shown in Sect. 4.3, we find values for β that are higher than −0.5, and this is accentuated from the analysis of the VEGA data (β ~ −0.3).

Apart from this issue regarding the analysis in Hα, we are able to describe well the disk density with the same physical parameters in both the Hα and Brγ lines: Σ₀ = 0.12 g cm⁻² and m = 3.0. As will be later discussed, this result found using HDUST is consistent with the ones presented in Sect. 4.3, showing a similar disk extension in these lines.

In Fig. 10, the intensity maps for each model are shown at the close-by continuum region and at different wavelength values in both the Hα and Brγ emission lines. The integrated intensity map (around each of these lines) is also presented. For a more realistic comparison, here we consider our best-fit kinematic model with a small flux contribution of 5% from the disk in the continuum nearby to Hα and a_c = 2 D_⋆. As shown in Table 2, these parameters were adopted as null in the kinematic analysis for the VEGA data, since we were not able to resolve the disk from our analysis of VEGA V² measurements in the continuum band (Sect. 3). Regarding the continuum region close to Brγ, the disk extension and flux contribution are given in Table 2 for the AMBER analysis.

The major difference between the intensity maps in Hα and Brγ is the disk flattening which is due to the different inclination angle derived from these two regions, i ~57° (Hα) and ~72° (Brγ), from the best models provided in Table 4. Moreover, as seen in the images, the stellar flattening is taken into account in the HDUST modeling, but not in the kinematic model (the star is modeled as a uniform disk). Apart from these departures, we see that our best-fit HDUST model presents a fairly similar distribution to the one computed with the kinematic code: a Gaussian distribution represents the circumstellar disk. This can be better noted considering the full integrated images around the emission lines.

7 Comparison to non-interferometric observables

In this section, we compare our best-fit models, found from the analysis of interferometric observables, to the observed spectral energy distribution (SED) and line profiles (Hα and Brγ) of o Aquarii. With respect to polarimetric data, it is discussed in Sect. 8.4.3 when addressing the disk stability.

Fig. 11 Comparison between the observed o Aquarii and model SEDs from the ultraviolet to the far-infrared region. Flux unit is in erg cm−2 s−1 Å−1 and wavelength is shown in logarithmic scale. IUE/SWP and IUE/LWP spectra are shown in black line and photometric data in black points. Top panel: purely photospheric models (color lines) with variation in the stellar radius (no inclusion of geometrical oblateness): R⋆ = 3.2 R⊙ (orchid), 4.0 R⊙ (red), and 4.4 R⊙ (green). Bottom panel: photospheric model with 4.0 R⊙ (red) and our best-fit HDUST model from fitting all the interferometric data (dashed blue line; Table 5). Note that the UBV-bands are better reproduced with R⋆ = 4.0–4.4 R⊙. Our best HDUST model reproduces the observed IR excess due to the circumstellar disk well.

7.1 Spectral energy distribution

In Fig. 11, we present the spectral energy distribution (SED) of o Aquarii from the ultraviolet (IUE/SWP and IUE/LWP spectra⁷) to the far-infrared region. References for the photometric data are given as follows: UBVJHK-bands (Anderson & Francis 2012), i-band (Henden et al. 2016), LM-bands (Bourges et al. 2017), and IRAS 12, 25, and 60 μm bands (Abrahamyan et al. 2015).

For the spectral region up to the V-band, we compare the data to the SEDs of purely photospheric atmosphere models with solar metallicity (Castelli & Kurucz 2004). In this region, the circumstellar disk flux level is much lower than the photospheric flux, thus allowing a proper probe of the stellar radius (e.g., Meilland et al. 2009). The surface gravity was fixed at log g = 4.0, this being the closest value in Castelli &Kurucz (2004) to logg = 3.9 that is given by our results of M_⋆ = 4.2 M_⊙ and R_⋆ = 4.0 R_⊙. The effective temperature was fixed at 13 000 K, following Cochetti et al. (2019). As in the previous sections, we consider the distance to be 144 pc, from the Gaia DR2 parallax.

These synthetic SEDs were calculated for three different stellar radius values, R_⋆ : 3.2 R_⊙ (Sigut et al. 2015), 4.0 R_⊙, and 4.4 R_⊙ (Cochetti et al. 2019). The value of 4.0 R_⊙ corresponds to the stellar radius determined from the fit to the VEGA V² data using a two-component model: 4.0 ± 0.3 R_⊙. The effect of interstellar medium extinction is not included in these models since it is negligible for o Aquarii. Assuming a total to selective extinction ratio of R_V = 3.1, Touhami et al. (2013) derived a color excess of E(B−V) = 0.015 ± 0.008 for this star from their fit to the SED. This means the observed flux is ~96% of the intrinsic one in the V-band (lower by ~0.02 dex). It is beyond the scope of this paper to estimate the extinction due to the circumstellar disk, however, from the comparison to purely photospheric models, we see in Fig. 11 that the effect of extinction (due to the interstellar and circumstellar matter) is conspicuously weak on the 0.220 μm bump.

From Fig. 11, we see that the UV and visible regions are better reproduced for a stellar radius of about 4.0–4.4 R_⊙, when compared to 3.2 R_⊙, adopted in Sigut et al. (2015), which corresponds to the expected polar radius for a B7 dwarf. We stress that the radius derived by Cochetti et al. (2019) is closer to our results from the fit to the VEGA V² data (Sect. 3). Their result of R_⋆ = 4.4 R_⊙ corresponds to a uniform disk diameter of θ ~ 0.28 mas (d = 144 pc). A better comparison to Cochetti et al. (2019) is hard since they do not provide error bars on R_⋆ from fitting the SED. Furthermore, they derived R_⋆ = 4.4 R_⊙ for o Aquarii using a distance of 134 pc from van Leeuwen (2007), rather than the value of 144 pc adopted here. From Fig. 11, this implies a larger discrepancy between the observed and synthetic SED for R_⋆ = 4.4 R_⊙, overestimating the observed flux.

We also compare the predicted SED of our best-fit HDUST model (Table 5) to the SED of the purely photospheric model with 4.0 R_⊙. Despite being able to reproduce the UBV-bands well, one sees that a purely photospheric model clearly underestimates the observed flux beyond the near-infrared due to the flux contribution from the circumstellar disk (e.g., Poeckert & Marlborough 1978; Waters 1986). From Fig. 11, it is evident that the SED is much better reproduced up to the far-infrared region when taking into account the IR excess from the gaseous circumstellar disk present in our best-fit HDUST model.

7.2 Hα and Brγ profiles

Our Hα spectra taken with the VEGA instrument (20 spectra, period from 2012 to 2016) are not analyzed in this work since they are saturated. This is a known effect seen in previous works on Be stars and correlated to the magnitude of the object. We stress that this instrumental saturation effect does not impact the visibilities and phases extracted from the fringes measured with VEGA (see, e.g., Delaa et al. 2011). To overcome this problem we used Hα line profiles from the BeSOS⁸ catalog (Arcos et al. 2018; Vanzi et al. 2012), obtained between 2012 and 2015, and thus covering a similar period to our VEGA observations. The typical spectral resolution of the BeSOS spectra is ~0.1 Å.

In Fig. 12, we compare the Hα and Brγ profiles from our best-fit models to observed profiles, namely, the mean Hα line profiles from BeSOS (7 profiles⁹) and the mean Brγ line profiles from our AMBER observations (8 profiles). The observed profiles in Fig. 12 were binned in wavelength in order to have a spectral resolution equal to one of the synthetic profiles from the kinematic and HDUST models: 1.3 Å (Hα) and 1.8 Å (Brγ). The mean EW in Hα from the BeSOS data is 19.1 Å. This is in agreement with the mean value of 19.9 Å found in Sigut et al. (2015), based on contemporaneous spectra, and adopted in our analysis with the kinematic code (Sect. 4.3).

First, we note that our best-fit kinematic and HDUST models provide a fairly reasonable match to the observed Hα and Brγ line profiles.The kinematic models correspond to our best-fits obtained from modeling the VEGA and AMBER differential data separately (Sect. 4). On the other hand, our best-fit HDUST model shown in Hα and Brγ is derived from the simultaneous fit to all our interferometric data (Table 5). Moreover, we stress the difficulty found by Sigut et al. (2015), using the radiative transfer code BEDISK, to reproduce the line wings and central absorption in the Hα profile of o Aquarii (see their Fig. 5).

However, it can be seen in Fig. 12 that both our best-fit kinematic and HDUST model are not able to properly reproduce, in particular, the wings of the Hα profile. On the other hand, the wings of the Brγ profile are fairly well reproduced by both of them, especially with HDUST.

Therefore, this inability to reproduce the wings of the Hα profile well is likely due to physical processes in the disk that are not taken into account in our models. It is known that the Hα profile wings of Be stars can be highly affected by non-coherent scattering, thus resulting in non-kinematic line-broadening in this transition (see, e.g., Hummel & Dachs 1992; Delaa et al. 2011). It is beyond the scope of this paper to quantify this possible effect in the Hα line of o Aquarii.

Fig. 12 Comparison between our best-fit kinematic models (dashed red; Table 2) and HDUST model (dashed blue, Table 5) in the Hα and Brγ line profiles. Mean observed line profiles of Hα (BeSOS) and Brγ (AMBER) are shown in black line. Our best-fit kinematic and HDUST models provide reasonable synthetic profiles to the observed ones in both Hα and Bγ.

8 Discussion

8.1 Disk extension in Hα and Brγ

In Sect. 4.3, we showed that the disk extension is similar in the Hα and Brγ lines. Interestingly, from previous studies, we could expect to find a larger disk extension in Hα than Brγ. For example, Meilland et al. (2011) found that δ Scorpii (B0.3IV), which was also observed with the VEGA and AMBER instruments, shows a circumstellar disk 1.65 times larger in Hα than in Brγ. Furthermore, Gies et al. (2007) derived the angular sizes of four Be stars (γ Cassiopeiae, ϕ Persei, ζ Tauri, and κ Draconis) in the K-band region using interferometric data from the CHARA/CLASSIC instrument. They showed that the disk of these stars was significantly larger (up to ~1.5–2.0 times) in the Hα line than in the K-band. However, Carciofi (2011) investigated theoretically, using the code HDUST, the formation loci of Hα and Brγ, and found them to be quite similar at least in the parameter space explored by the authors (see their Fig. 1). Moreover, Stee & Bittar (2001), using the code SIMECA, found that Be star disks can be larger (up to two times) in Brγ than in Hα.

For a quantitative comparison of the disk extension in Hα and Brγ, we fitted simple Gaussian distributions to the intensity map of our best-fit HDUST model for all the values of inclination angle in BeAtlas. In order to remove the contribution from the star and disk continuum, we removed the image from the continuum before performing the fit and we hide the central part of the image which is affected by the stellar contribution.

In Fig. 13, we show the major-axis FWHM from our fit as a function of the inclination angle for the Hα and Brγ lines. First, one sees that the disk size-extension (major-axis FWHM) varies differently in the Hα and Brγ lines as a function of the inclination angle. The disk extension increases in Brγ with the inclination angle. On the other hand, it decreases significantly in Hα up to i ~ 56° and increases after this value. One sees that the ratio between the extension in these lines decreases from about 1.50 at zero inclination to about 1.05 at 63.5°. Furthermore, we note that the disk extensions in these lines are very close to each other for i ~ 56° (Hα) and i ~ 72° (Brγ): major-axis FWHM ~ 2.45 mas. Considering d = 144 pc, the disk size is ~10 D_⋆ (close to our findings from the kinematic modeling).

Therefore, from this simple analysis using HDUST models, we verify our findings using the kinematic model: a similar circumstellar disk extension in Hα and Brγ. This arises since the (equivalent) Gaussian disk to our best-fit HDUST model presents quite different changes on its extension in these lines as a function of the inclination angle. Based on that, we can also explain the difference between δ Scorpii and o Aquarii. The former is seen under a low inclination angle (~30°) and exhibits a high ratio between the Hα and Brγ disk sizes. The latter is seen under a higher inclination angle and shows similar disk sizes in both lines. On the other hand, as discussed above, ϕ Persei and ζ Tau show larger disks in Hα than in the K-band and these stars are seen close to edge-on with i = 78° (Mourard et al. 2015) and 85° (Carciofi et al. 2009), respectively. Thus, this similarity in the disk extensions, found for quite different values of inclination angle, could indicate a more complex physical structure of the circumstellar disk than the one assumed by our best-fit HDUST model (based on a vertically isothermal disk).

Fig. 13 Major-axis FWHM of Gaussian distribution (fitted from our best-fit HDUST model) as a function of the HDUST inclination angle. All the other HDUST parameters are fixed. Blue points correspond to the fit in Hα and red points in Brγ. The vertical dashed lines mark our values for inclination angle derived from the HDUST analysis, fitting the data in Hα (blue) and Brγ (red). Note that the equivalent Gaussian fits show a similar extension (2.45 mas, marked in horizontal dashed line) for these values of i.

8.2 Inclination angle and vertical disk structure

From our Hα and Brγ differential data analysis, using the kinematic model, we achieved good precision in the determination of the stellar inclination angle: i ~ 61.2 ± 1.8° (VEGA) and i = 75.9 ± 0.4° (AMBER). Nevertheless, there is a clear discrepancy between the inclination angle found from fitting the VEGA and AMBER datasets. The value determined from VEGA is about 15° lower than the one found in the analysis of the AMBER data. We can show that this issue does not stem from an intrinsic limitation of the kinematic code (2D model) for Be stars seen under high inclination angle (i ≳ 60°). Indeed, by using a sophisticated 3D radiative transfer model (HDUST), not subjected to such a limitation, we verified the same discrepancy on i from the fit to these datasets separately (see, again, in Fig. 9, the trend of ${\chi }_{\text{r}}^2$, as a function of i).

It may be argued that the difference found in inclination angle is hiding a difference in the disk thickness in these lines. Assuming a non-geometrically thin disk, for an ellipse with major and minor axes denoted, respectively, by a and b, the ratio between a and b, the circumstellar disk flattening, is given by (see, e.g., Meilland et al. 2007b): $$\frac{a}{b}=\frac{1}{\cos i+2\sin \frac{\Theta }{4}\sin \left(i-\frac{\Theta }{4}\right)},$$(10)

where i is the stellar inclination angle and Θ the disk opening angle. Since the i derived from Hα using our physical models is much lower than from Brγ (a reliable value when compared to other results in the literature), this would imply a disk thicker (higher opening angle Θ) in Hα than in Brγ. Considering the values described above, the disk opening angle in Hα would be Θ ~ 37° larger in Hα than in Brγ (assuming a geometrically thin disk in Brγ). Such a high value of opening angle is far beyond what is measured and expected by the VDD model, typically less than ~10° (cf. Rivinius et al. 2013). This might indicate the necessity of more complex physical assumptions in the physical properties of our disk model.

Since the code HDUST provides a pure hydrogen modeling for the photosphere plus disk regions, this disagreement between the VEGA (Hα) and AMBER (Brγ) analyses in the determination of i could be due to an opacity effect. It is well-known that the inclusion of heavy elements can impact the density and temperature stratifications in the circumstellar disk of Be stars by shielding emission from the star. (see, e.g., Sigut & Jones 2007). Furthermore, we stress that our best-fit HDUST model is a parametric model (based on a vertically isothermal structure). Departures from vertically isothermal disks are well-known in the literature. For example, using the radiative transfer code BEDISK, Sigut et al. (2009) verified that isothermal and self-consistent hydrostatic models can present large differences regarding the temperature stratification in the disk of Be stars. Using HDUST, Carciofi & Bjorkman (2008) also found that non-isothermal effects can be significant for denser Be star disks. Thus, further investigation is needed concerning this effect on the determination of i for o Aquarii, but that is beyond the scope of this paper.

Finally, another possibility to explain the difference in apparent inclination angle found in our modeling could be a non-negligible contribution of a polar wind. Clues of the presence of polar wind, or at least of circumstellar material in the polar regions, have been found by Kervella & Domiciano de Souza (2006) and Meilland et al. (2007a). In our models, we assume that all the circumstellar material is in the thin equatorial disk. If a non negligible fraction of the material is located near the poles, although we would expect it to be quite diluted and optically thin (at least in the continuum), it might affect the line emission with a different magnitude in Hα and in Brγ. If one assumes that the hydrogen level populations favor Hα emission over Brγ, the polar contribution of Hα would be higher, and the environment might look less flattened in this line than in Brγ.

8.3 Stellar and disk rotation

In Sect. 5, our results are presented in terms of the stellar obla- teness R_eq /R_p (denoted by f in Eq. (11)). First, we give the relation between the oblateness and the angular Ω∕Ω_crit and linear v_rot∕v_crit rotational rates as follows: $$\frac{\Omega }{{\Omega }_{\text{crit}}}=\frac{v_{\text{rot}}}{v_{\text{crit}}}\frac{R_{\text{eq},\text{crit}}}{R_{\text{eq}}}={\left(\frac{3}{2}\right)}^{3/2}{\left[\frac{2(f-1)}{f^3}\right]}^{1/2},$$(11)

where R_eq,crit and R_eq (in units of polar radius) are, respectively, the stellar equatorial radius in the case of critical velocity and the actual one (see, e.g., Frémat et al. 2005; Ekström et al. 2008).

Considering only the uncertainties on v_rot (Table 2), with the critical velocity v_crit fixed to 391 km s⁻¹ (Frémat et al. 2005), we obtain a linear rotational rate of v_rot /v_crit = 0.83 ± 0.02 (v_rot = 325 ± 6 km s⁻¹, VEGA) and 0.775 ± 0.005 (v_rot = 303 ± 2 km s⁻¹, AMBER). From the HDUST analysis, we find v_rot /v_crit = 0.96 (VEGA and AMBER) from our best-fit model (no error bars). This difference between the kinematic and HDUST analysis can be explained since the β exponent (velocity law in the disk) is fixed in the HDUST analysis (Keplerian disk, β = −0.5), while it is a free parameter in the kinematic model. We derived values for β from the kinematic analysis that are significantly higher (more positive) than −0.5 (see Table 2).

Apart from these differences, our analysis is consistent with a high rotational rate for o Aquarii, showing v_rot /v_crit from ~0.8 up to 1.0, depending on the particular analysis considered. The BeAtlas fits to the VEGA and AMBER differential data are significantly worsened (Fig. 9), when considering R_eq /R_p = 1.20–1.30 (Ω∕Ω_crit = 0.88–0.96). Thus, our HDUST analysis indicates that o Aquarii rotates faster than Ω∕Ω_crit = 0.96, disfavouring the lower range of Ω∕Ω_crit between 0.86 and 0.93 that is derived by Cochetti et al. (2019).

In Sect. 4.3, we found a strong correlation between the velocity at the base of the disk and the β exponent of the rotation law β. The inferred degeneracy, stronger in the case of the VEGA data, which have a lower spectral-resolution with higher uncertainties, prevents us from independently constraining these two parameters with our kinematic model when fittingonly our spectro-interferometric data. However, the addition of an external constraint, the measured v sin i, removed this degeneracy, allowing us to derive a more accurate value of β in comparison to the other MCMC fitting tests (Appendix B).

From our MCMC fit to the AMBER dataset, with a preset v sin i, we derived a β of −0.426 ± 0.003. Thus, the disk appears to be rotating in a nearly Keplerian fashion. Despite this very low error on β, note that the error bars on β change with respect to the presented MCMC tests, up to ± 0.008 (see Fig. B.2). On the other hand, the value derived from the fit of the VEGA data is about 0.1 higher than from AMBER. We stress that this discrepancy cannot be explained by a radial dependent rotational law because both lines roughly stem from the same region in the disk (similar disk extensions in these lines).

Moreover, this apparent higher value of β in Hα was also theorigin of some biases that we found when modeling the VEGA differential data alone using HDUST. Without fixing M_⋆, the VEGA analysis with HDUST favours unrealistically high values of stellar mass up to ~11 M_⊙. This happened due to the fact the higher mass models also correspond to higher rotational velocity at the base of the disk. As the value of β is fixed to −0.5 in the BeAtlas grid of models, this was the only way to increase the rotational velocity in the disk. In Fig. 14, our best-fit kinematic and HDUST models are compared to a higher mass HDUST model for one VEGA differential measurement and the Hα profile. When compared to our best HDUST model (mass fixed to 4.2 M_⊙), the visibilitydrop and the Hα profile are better reproduced with HDUST models with higher value of mass, but also considering a larger value of Σ₀ (0.28 g cm⁻²). This happens since the stellar radius is also increased for a higher mass model and the flux contribution from the star is larger in the line. In this case, our BeAtlas model is able to produce more similar synthetic Hα visibility and profile to the ones from our best-fit kinematic model.

One possible explanation for the discrepancy between the value of β determined from Hα and Brγ could be the higher effects of non-kinematic broadening on Hα. This is already evidenced by the larger wings, in terms of Doppler shift, for this emission line. Such effects are known to be due to non-coherent scattering in the circumstellar environment, as explained, for example, in Auer & Mihalas (1968). Global effects on interferometric data were discussed by Stee et al. (2012) in the case of the Be star γ Cassiopeiae observed with VEGA. These authors used a similar kinematic model, but with two additional parameters to quantify the non-coherent scattering and found that about half the flux in the line was affected by such an effect. Never- theless, the possible bias on the measurement of β in a line strongly affected by such non-kinematic broadening should be investigated further.

We also note that a possible close companion could influence v_rot, as well as the disk structure, as previously mentioned. However, the presence of a close companion with a detectable influence on the measured parameters seems excluded from the observed calibrated V², and in particularfrom the spectro-interferometric differential observables, which both show signatures well reproduced by a symmetric rotating disk.

Fig. 14 Bias effect of the disk velocity on Hα modeling. One VEGA measurement and observed Hα profile (BeSOS, as in Fig. 12) are shown in black lines. Our best-fit kinematic (dashed red) and HDUST (dashed blue) models are shown in Hα visibility and line profile. They are compared to HDUST models with a higher mass of 10.8 M⊙ with: Σ0 = 0.12 g cm−2 (dashed orchid) and Σ0 = 0.28 g cm−2 (dashed green). See text for discussion.

Fig. 15 Top panel: 70 observed Hα profiles of o Aquarii from the BeSS database, covering about 17 yr of observations (2001–2018). Bottom panel: Hα equivalent width as a function of the observation time (modified Julian date). Civil dates are indicated for a part of the measured EW. The mean EW (solid) within the standard deviation (dashed lines) is marked in black. Local regression fit of EW, as a function of time, is shown as a solid red line. The time interval covered by our interferometric observations (VEGA) is indicated with dashed red lines. See text for discussion.

8.4 Disk variability: a multi-technique analysis

8.4.1 Spectroscopy

The Be star o Aquarii is known to possess a stable Hα line profile for up to several years. For example, Sigut et al. (2015) verified that the EW in the Hα line is stable (within about 5%) up to about nine years (from 2005 to 2014).

To go further in the analysis of the disk stability, we analyzed 70 Hα line profiles, spanning from 2001 to 2018, from the BeSS database¹⁰ (Neiner et al. 2011). Since these observations are performed with several instruments, the line profiles shown here are interpolated to have spectral resolution of 0.5 Å (lowest resolution in the dataset). From these observations, we calculated the equivalent width (EW) in the Hα line. In Fig. 15, we show the analyzed Hα spectra together to the temporal evolution of the Hα EW.

We found that the disk is fairly stable over this 17-yr time span with a mean value of $\overline{\text{EW}}$ = 18.1 ± 1.2 Å. This valueagrees well to older results in the literature. Slettebak & Reynolds (1978) measured Hα EW = 18.80 ± 0.11 Å in 1975 and 18.58 ± 0.21 Å in 1976. From the Hα profile observed in 1981, Andrillat (1983) measured EW = 17.2 Å. Thus, this supports an even longer global disk stability up to at least 40 yr. However, a slight increasing trend in EW is seen between 2001 and 2012. This could suggest an augmentation in the disk density of o Aquarii in this period. Considering the period of our interferometric observations (from 2012 to 2016), it is hard to observe any trend of Hα EW as a function of time.

Fig. 16 VEGA (top panels) and AMBER (bottom panels) differential visibilities extracted from observations at different epochs (black line). VEGA measurements span four years and the AMBER ones span three years. The observation date and the baseline length (projected onto the sky) are indicated in the top of each panel. Our best-fit kinematic models derived from the fit, in a separate way, to each dataset (VEGA and AMBER) are shown in dashed red line. Note that our best-fit kinematic models match well to the differential visibilities obtained at different epochs.

8.4.2 Interferometry

These results are consistent with our ability to model, with the same model parameters, simultaneously all our VEGA and AMBER data regardless of the epoch. In Fig. 16, we present a temporal evaluation of our spectro-interferometric data in the Hα and Brγ lines. Since the drop in visibility is expected to change due to possible variations in the disk extension, we only show here the differential visibilities from the VEGA and AMBER observations. These measurements are chosen to cover the whole period of our observations from 2011 to 2016. For a more robust comparison, we chose measurements obtained with different baseline lengths (projected onto the sky), and, thus, covering different levels of spatial resolution.

One sees that, regardless of the period of time of the observations, our final kinematic models provide a very reasonable match to both the VEGA and AMBER data. Thus, considering our interferometric data, we are not able to detect any conspicuous variation of the circumstellar disk extension within a period up to five years (from 2011 to 2016). This is in agreement with previous interferometric studies of o Aquarii by Sigut et al. (2015). Besides that, this analysis supports our approach of fitting each one of the interferometric datasets (VEGA and AMBER) without imposing any discrimination based on the observation time.

Fig. 17 Polarimetric quantities of o Aquarii, as a function of the observation time, spanning about six years. Top panel: observed V-band polarization (44 measurements). Bottom panel: ratio between the observed B- and R-bands polarization (40 measurements). The mean values of these quantities are shown in dashed line. See text for discussion.

8.4.3 Polarimetry

Additional multi-epoch polarimetric data also support our findings of a stable disk for o Aquarii, close to the steady-state regime. In Fig 17, we show the temporal evolution of broad-band linear polarimetry in the V-band (P_V) of o Aquarii, as well as the ratio between the B- and R-bands polarization (P_B∕P_R).

These data were obtained over 43 nights, from June 2010 to August 2016, with the IAGPOL polarimeter (Magalhães et al. 1996), mounted on the 0.6 m Boller & Chivens telescope at Observatório do Pico dos Dias (OPD/LNA). This polarimeter is composedby a rotating half-wave retarder and a Savart Plate used as analyser to provide the modulation of the light polarization, and then the polarimetric quantities. Details of data reduction are found in Magalhães et al. (1984) and Bednarski (2016).

From Fig. 17, the mean value of the observed V-band polarization is $\overline{P_V}$ = 0.48 ± 0.03%. This value, derived from the mean and standard deviation of the Stokes Q and U parameters, is compatible to the one determined by Yudin (2001), namely, 0.52 ± 0.05%. Since the observations from Yudin (2001) predate our OPD/LNA observations by more than a decade, we conclude that the polarization values of o Aquarii remained very constant for over 20 yr.

In order to determine the intrinsic value of polarization, the interstellar contribution to the observed values quoted above must be removed. For that, we observed four main sequence stars, in the BVRI-bands, which are angularly close to o Aquarii. A MCMC method was implemented to process the four BVRI data of each field star, generating a sample of the likelihood function in terms of the interstellar Serkowski parameters P_max and λ_max (Serkowski et al. 1975; Wilking et al. 1982). The best estimates for these parameters are shown in Table D.1.

There is a good agreement among the PA values of the field stars. Moreover, by using Gaia DR2 distances, we found that P_max increases linearly along the line of sight of o Aquarii (seeFig. D.1). In this case, it suggests that the alignment of the grains at the interstellar medium is nearly homogeneous (e.g., McLean & Clarke 1979). Thus, from a simple linear fit to P_max vs. distance for the field fields, we determined P_max for o Aquarii. The derived interstellar polarization parameters for o Aquarii are shown in Table 6, which are in reasonable agreement with the ones reported in Yudin (2001) of P_max = 0.20% and PA_IS = 125° (no error bars).

Taking into account our results for the interstellar polarization components, we found for the intrinsic V-band polarization and position angle $P_V^{\text{int}}$ = 0.49 ± 0.03% and PA ^int = 2.5 ± 2.7°, respectively. Yudin (2001) determined $P_V^{\text{int}}$ = 0.60% with PA^int = 6.0°, which is close to our PA^int value. Moreover, both estimates for PA^int are consistent with our determination for the disk major-axis position angle (~110°), being almost perpendicular to the polarization vector, as expected.

Furthermore, our best-fit HDUST model (Table 5) predicts a polarization degree of 0.41% in the V-band. This agrees well with our measurement for the average intrinsic polarization of the OPD/LNA data. Therefore, besides the independent checks provided by the SED and spectroscopic data (Sect. 7), our polarimetric data also support our physical model for o Aquarii, which was derived purely from the fit to interferometric data (as discussed in Sect. 5.3).

Lastly, Fig. 17 shows that both the polarization degree in the V-band and the ratio between the B- and R-bands are almost constant in time, showing a small scatter around the mean value. In particular, this latter quantity is related to the density scale at the inner portion of the disk (Haubois et al. 2014). From the theoretical investigation of Panoglou et al. (2019), the variation on the polarization degree in the V-band (ΔP_V) can reach up to about 0.1% due to asymmetries in the disk density structure, caused by a binary companion. Moreover, Haubois et al. (2014) predicted ΔP_V of up to 2% due to temporal changes in the mass decretion rate. The standard deviation of our P_V distribution (approximately Gaussian), namely, ~0.03%, is quite a bit lower than the above values. It is well explained in terms of the precision of our polarimetric data, as the typical error bar on P_V is ~0.01–0.02% (Fig. 17).

8.4.4 A stable disk

Besides the analysis of the Hα EW and broad-band polarimetric quantities, our modeling with the code HDUST indicates that the disk must be close to the steady-state regime: having a radial density law exponent of 3.0 (e.g. Haubois et al. 2012; Vieira et al. 2017). Other studies of o Aquarii are in fair agreement to our findings from HDUST. Using the radiative transfer code BEDISK, Silaj et al. (2010) derived m = 3.5 from the fit to the Hα profile, while Sigut et al. (2015) found m = 2.7 as a representative value from the analysis of all the different observables.

Previous and ongoing studies of Be stars with stable disks found similar results to ours. For example, Klement et al. (2015) found m = 2.9 for the late-type Be star β Canis Minoris (B8Ve). Correia Mota (2019) derived m = ${2.44}_{-0.16}^{+0.27}$ for α Arae (B2Vne). The B9Ve star α Columbae shows m = ${2.54}_{-0.13}^{+0.06}$ (A. Rubio, priv. comm.). Thus, the radial density exponent is consistently equal or somewhat less than 3.0 for these Be stars with stable disks.Also, from analysing the temporal variation of the disk density, Vieira et al. (2017) identified a slightly extended range of m (between ~3.0 and ~3.5) for the steady-state regime, in comparison to the canonical value of 3.5. As pointed out by these authors, this canonical value is based on simplifications of the standard theory, which assumes, for example, vertically isothermal disks and isolated systems (single stars). One possibility to explain the measured m lower than 3.5 could thus rely on non-isothermal effects in the disk structure (see, e.g., Carciofi & Bjorkman 2008).

Finally, we note that such long-term stability of o Aquarii’s disk is consistent with other results in the literature: late-type Be stars are more likely to have more stable disks than earlier Be stars (e.g., Vieira et al. 2017; Labadie-Bartz et al. 2018; Rímulo et al. 2018). As discussed in Sect. 8.3, the stellar rotation seems to be very close (~96%) to the critical value (391 ± 27 km s⁻¹ from Frémat et al. 2005), in particular regarding the HDUST analysis: v_rot = 368 km s⁻¹ (Table 5). This is consistent with the results from Cranmer (2005): Be stars with lower effective temperature T_eff ≲ 21 000 K – that is, later spectral types such as our target – are more likely to have a rotation rate close to one than the earlier Be stars. Thus, one possibility to explain such a long-term stability of the disk of o Aquarii could rely on its fast rotation, ensuring in this case a nearly constant mass-injection rate into the disk.

9 Conclusions

We analyzed VEGA V², as well as VEGA and AMBER differential visibility and phase of the Be-shell star o Aquarii. To date, the spectro-interferometric dataset analyzed in this paper is the largest for a Be star, considering quasi-contemporaneous observations in both the Hα (VEGA) and Brγ (AMBER) lines.

For the first time, we measured o Aquarii’s stellar radius (R_⋆ = 4.0 ± 0.3 R_⊙) and determined the disk extension in the Hα and Brγ lines as, respectively, 10.5 ± 0.3 D_⋆ and 11.5 ± 0.1 D_⋆. Using radiative transfer models computed with the code HDUST, we explained the quasi-identical extension of the emission in these lines by an opacity effect found for disks seen under a high inclination angle.

We showed that the inclination angle derived from Hα is about 15° lower than the one determined in Brγ, when analysing each line separately with HDUST. More complex physical models, for example, with non-isothermal vertical scaling of the disk or the addition of heavier elements, could resolve this issue and should be investigated in the future.

Our simple kinematic model highlighted the high correlation between the rotational velocity at the base of the disk and the rotational law exponent β. Assuming external constraints, such as v sin i, we managed to constrain this parameter and showed that the disk rotation is nearly Keplerian (β ~ 0.43) from the analysis in the Brγ emission line. As for the inclination angle, the determination of β, using the Hα line (β ~ 0.30), seemsto be significantly biased. Other studies also verified such a large deviation from the Keplerian rotation for Be stars when analysing interferometric quantities measured in Hα (see, e.g., Delaa et al. 2011). One possible explanation would be the higher effect of non-coherent scattering on the Hα line formation than on Brγ.

Despite being derived purely from the fit to interferometric data, our best-fit HDUST model provides a very reasonable match to non-interferometric observables of o Aquarii: the observedSED, Hα and Brγ line profiles,and polarimetric quantities. Thus, this cross-check provides an independent validation of our best-fit physical model. We found using HDUST a satisfying common physical description for the circumstellar disk in both Hα and Brγ: a base disk surface density Σ₀ = 0.12 g cm⁻² (ρ₀ = 5.0 × 10⁻¹² g cm⁻³) and a radial density law exponent m = 3.0, that is, close to the steady-state regime according to the VDD model (m = 3.5). This result agrees with recent studies of other Be stars with stable disks, and may indicate the necessity to revise m = 3.5 (steady-state standing for single stars with vertically isothermal disks) that is predicted by the VDD theory. Otherwise, this could indicate non-isothermal effects on the disk vertical structure of o Aquarii. The long-term stability of the o Aquarii’s disk is verified by our analysis of a large sample of Hα profiles and polarimetric data, spanning about 20 and six years, respectively. Combined with older results in the literature, a longer global disk stability is suggested for up to at least 40 yr.

The stellar rotation seems to be very close (~96%) to the critical value (391 km s⁻¹), in particular accordingly to our HDUST analysis: v_rot = 368 km s⁻¹ from the best-fit HDUST model with fixed M_⋆ = 4.2 M_⊙ (cf., Sects. 5.2 and 5.3). One possibility to explain such a long-term stability in the disk of o Aquarii could rely on its own high stellar rotation, being, in this case, a main source for the mass injection from the stellar surface to the disk. Thus, apart from the mass decretion due to other possible mechanisms in Be stars, this would provide a constant rate of mass injection. In short, our results on the stellar rotation and on the disk stability are consistent with the literature results showing that late-type Be stars are more likely to be fast rotators and have stable disks (see Sect. 8.4.4).

Finally, to further investigate these issues, our multi-wavelength and multi-emission line modeling approach must be performed on a larger sample of Be stars with disks of different densities and seen under different inclination angles. The implementation of a MCMC model fitting procedure with the kinematic model, and the use of our grid of HDUST models (BeAtlas), are very promising for the spectro-interferometric analysis of a large survey of Be stars, providing robust model parameters and associated uncertainties. A future project will attempt this task on a few dozen objects observed with VEGA and AMBER.

Acknowledgements

We thank the anonymous referee for helping to improve this paper. E.S.G.de.A. thanks OCA and the “Ville de Nice” (Nice, France) for the financial support to this work through the “Bourse Doctorale Olivier Chesneau” during the period of 2016–2019. E.S.G.de.A. acknowledges A. Rubio for relevant information about her work on α Columbae. R.L. has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement n. 664931. D.M.F acknowledges FAPESP (grant 2016/16844-1). A.C.C. acknowledges support from CNPq (grant 307594/2015-7). This work was supported by the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU co-funded by CEA and CNES. This work 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 work is based upon observationsobtained with the Georgia State University Center for High Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA Array is supported by the National Science Foundation under Grants No. AST-1211929 and AST-1411654. This work used BeSOS Catalogue, operated by the Instituto de Física y Astronomía, Universidad de Valparaíso, Chile: http://besos.ifa.uv.cl and funded by Fondecyt iniciación N 11130702. Based on observations collected at the European Southern Observatory under ESO programmes 087.D-0311 and 094.D-0140. This work has made use of the BeSS database, operated at LESIA, Observatoire de Meudon, France: http://basebe.obspm.fr. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This research has made use of the Jean-Marie Mariotti Center (JMMC) services LITpro, SearchCal, and AMHRA, co-developped by CRAL, IPAG and Lagrange. This work has made use of the SIMBAD and VizieR databases, operated at CDS, Strasbourg, France.

Appendix A Observational logs

Appendix B MCMC fitting tests: fits to the VEGA and AMBER data with the kinematic model

Fig. B.1 As in Fig. 4, but for the other MCMC fitting tests (test i in the left and test ii in the right) to fit the VEGA differential data.

Fig. B.2 As in Fig. 4, but for the other MCMC fitting tests (test i in the left and test ii in the right) to fit the AMBER differential data.

Appendix C Best-fit kinematic and HDUST models: AMBER

Fig. C.1 Comparison between our best-fit kinematic (dashed red; Table 2) and HDUST (dashed blue; Table 5) models to all the AMBER measurements (black).

Appendix D Interstellar polarization

Fig. D.1 Fitted Pmax for the field stars (Table D.1, open triangles) as a function of the Gaia DR2 distance. From the linear fit to Pmax vs. d for the field stars (dotted line), we determined Pmax for o Aquarii (redcross).

References

All Tables

All Figures

	Fig. 1 uv plan coverage obtained around Hα (0.656 μm) with CHARA/VEGA (green) and Brγ (2.166 μm) with VLTI/AMBER (red).
In the text

	Fig. 2 VEGA V2 measurements in the close-by continuum band (top) and in the Hα band (bottom) are shown in red points. Our best-fit models consisting of one (solid line) and two (dashed line) uniform disks are overplotted in blue. See Table 1 and text for discussion.
In the text

Fig. 3 Top panels: uniform disk diameter derived from each individual VEGA V2 measurements (continuum band in the left and Hα band in the right) plotted as a function of the baseline position angle (PA). The red dotted line represents the best-fit diameter from modeling all the data in each band (θ = 0.28 mas in the continuum and θ = 0.36 mas in the Hα band). Bottom panels: corresponding normalized residuals.

In the text

Fig. 4 Histogram distributions and two-by-two correlations (after the burn-in phase) for the free parameters of our best-fit kinematic models using MCMC for the VEGA (left panel) and AMBER (right panel) differential data. The median values are shown in solid red lines and the first and third quartiles in dashed red lines. The median and the first and third quartiles estimated for the parameters of our best-fit models (VEGA and AMBER) are presented inTable 2. In the correlation plots, darker points correspond to models with lower values of χ2. See text for discussion.

In the text

	Fig. 5 Comparison between our best-fit kinematic models (dashed red; Table 2) and two different VEGA (top panels) and AMBER (bottom panels) measurements (black line). Our best-fit HDUST model is also shown (dashed blue; Table 5; discussion in Sect. 6). δλ of the kinematic model and AMBER data is increased to 1.8 Å in order to compare them to the HDUST model (δλ fixed to 1.8 Å).
In the text

Fig. 6 ${\chi }_{\text{r}}^2$ maps of 40 000 kinematic models as a function of vrot and β from the fit to VEGA (top panel) and AMBER (bottom panel) differential data. Only these two parameters were varied in a regular step in the intervals shown here. The other parameters are fixed (Table 2). Our results found from the MCMC analysis for vrot and β are indicated with red crosses. In order to highlight the correlation between β and vrot, the gray region corresponds to an arbitrary number of models, encompassing about the 5000 best models in both cases. The value of β = −0.5 (Keplerian disk) and our determination for vrot are marked in dashed black line. Note the strong correlation between the stellar rotational velocity and the disk velocity law exponent in both the cases. Also, note that a Keplerian disk is found from modeling the AMBER data, but not from VEGA.

In the text

Fig. 7 Analysis (i): lowest value of reduced χ2 (${\chi }_{\text{r}}^2$) for each value of Req∕Rp (left panel) and stellar mass (middle panel) from the HDUST fit to the VEGA V2 data (642.5 nm band). Local regression fits to ${\chi }_{\text{r}}^2$, as a function of the parameter values, are shown in red line. In the right panel, the predicted visibility from our best-fit HDUST model (red points; Table 5) is compared to the VEGA V2 measurements in the continuum band (black points). The predicted visibility from the HDUST model with the highest mass in the BeAtlas grid (14.6 M⊙, highest ${\chi }_{\text{r}}^2$ in the middle panel) is shown in blue points.

In the text

	Fig. 8 Lowest value of ${\chi }_{\text{r}}^2$ for each value of disk major-axis position angle from the HDUST fit to the VEGA (top panel, analysis (ii)) and AMBER (bottom panel, analysis (iii)) differential visibility and phase. Local regression fits of ${\chi }_{\text{r}}^2$ as a function of the disk PA are shown as a red line.
In the text

Fig. 9 Lowest value of ${\chi }_{\text{r}}^2$ for each value of stellar inclination angle, oblateness, base disk surface density, and disk density law exponent from the HDUST fit to the: VEGA differential data (top, analysis (ii)), AMBER differential data (middle, analysis (iii)), and all the interferometric data considered in this section (bottom panel, analysis (iv)). The stellar mass is fixed to 4.2 M⊙ and disk PA to 110°. Local regression fits to ${\chi }_{\text{r}}^2$, as a function of the parameter values, are shown as a red line. The mean parameter values for the sets of best models (Table 4) are marked in dashed black line. Our best-fit BeAtlas model to fit all the interferometric data is shown in Table 5. See text for discussion.

In the text

	Fig. 10 Intensity maps of our best-fit HDUST and kinematic models at different wavelengths around the Hα line (first two rows) and the Brγ line (last two rows). Flux/pixel is in arbitrary units with the same scale in Hα and Brγ. The image integrated in wavelength around each of these lines (Δλ = 2.7 nm around Hα and 3.9 nm around Brγ) are shown in the second column.
In the text

Fig. 11 Comparison between the observed o Aquarii and model SEDs from the ultraviolet to the far-infrared region. Flux unit is in erg cm−2 s−1 Å−1 and wavelength is shown in logarithmic scale. IUE/SWP and IUE/LWP spectra are shown in black line and photometric data in black points. Top panel: purely photospheric models (color lines) with variation in the stellar radius (no inclusion of geometrical oblateness): R⋆ = 3.2 R⊙ (orchid), 4.0 R⊙ (red), and 4.4 R⊙ (green). Bottom panel: photospheric model with 4.0 R⊙ (red) and our best-fit HDUST model from fitting all the interferometric data (dashed blue line; Table 5). Note that the UBV-bands are better reproduced with R⋆ = 4.0–4.4 R⊙. Our best HDUST model reproduces the observed IR excess due to the circumstellar disk well.

In the text

	Fig. 12 Comparison between our best-fit kinematic models (dashed red; Table 2) and HDUST model (dashed blue, Table 5) in the Hα and Brγ line profiles. Mean observed line profiles of Hα (BeSOS) and Brγ (AMBER) are shown in black line. Our best-fit kinematic and HDUST models provide reasonable synthetic profiles to the observed ones in both Hα and Bγ.
In the text

Fig. 13 Major-axis FWHM of Gaussian distribution (fitted from our best-fit HDUST model) as a function of the HDUST inclination angle. All the other HDUST parameters are fixed. Blue points correspond to the fit in Hα and red points in Brγ. The vertical dashed lines mark our values for inclination angle derived from the HDUST analysis, fitting the data in Hα (blue) and Brγ (red). Note that the equivalent Gaussian fits show a similar extension (2.45 mas, marked in horizontal dashed line) for these values of i.

In the text

Fig. 14 Bias effect of the disk velocity on Hα modeling. One VEGA measurement and observed Hα profile (BeSOS, as in Fig. 12) are shown in black lines. Our best-fit kinematic (dashed red) and HDUST (dashed blue) models are shown in Hα visibility and line profile. They are compared to HDUST models with a higher mass of 10.8 M⊙ with: Σ0 = 0.12 g cm−2 (dashed orchid) and Σ0 = 0.28 g cm−2 (dashed green). See text for discussion.

In the text

Fig. 15 Top panel: 70 observed Hα profiles of o Aquarii from the BeSS database, covering about 17 yr of observations (2001–2018). Bottom panel: Hα equivalent width as a function of the observation time (modified Julian date). Civil dates are indicated for a part of the measured EW. The mean EW (solid) within the standard deviation (dashed lines) is marked in black. Local regression fit of EW, as a function of time, is shown as a solid red line. The time interval covered by our interferometric observations (VEGA) is indicated with dashed red lines. See text for discussion.

In the text

Fig. 16 VEGA (top panels) and AMBER (bottom panels) differential visibilities extracted from observations at different epochs (black line). VEGA measurements span four years and the AMBER ones span three years. The observation date and the baseline length (projected onto the sky) are indicated in the top of each panel. Our best-fit kinematic models derived from the fit, in a separate way, to each dataset (VEGA and AMBER) are shown in dashed red line. Note that our best-fit kinematic models match well to the differential visibilities obtained at different epochs.

In the text

	Fig. 17 Polarimetric quantities of o Aquarii, as a function of the observation time, spanning about six years. Top panel: observed V-band polarization (44 measurements). Bottom panel: ratio between the observed B- and R-bands polarization (40 measurements). The mean values of these quantities are shown in dashed line. See text for discussion.
In the text

	Fig. B.1 As in Fig. 4, but for the other MCMC fitting tests (test i in the left and test ii in the right) to fit the VEGA differential data.
In the text

	Fig. B.2 As in Fig. 4, but for the other MCMC fitting tests (test i in the left and test ii in the right) to fit the AMBER differential data.
In the text

	Fig. C.1 Comparison between our best-fit kinematic (dashed red; Table 2) and HDUST (dashed blue; Table 5) models to all the AMBER measurements (black).
In the text

	Fig. D.1 Fitted Pmax for the field stars (Table D.1, open triangles) as a function of the Gaia DR2 distance. From the linear fit to Pmax vs. d for the field stars (dotted line), we determined Pmax for o Aquarii (redcross).
In the text