© The Authors 2025

1 Introduction

The star π Aquarii (HD 212571, HIP 110672) is a bright rapidly rotating classical Be star (B1III-IVe) presenting mid- and longterm variability. Notably, the star has been observed for more than 100 years. McLaughlin (1962) reported variations over 50 years in the line profile of Hα and Hγ, where more of the profiles were observed in double-peaked line emission, probably indicating the active phase of the disk of π Aquarii. He also noted a significant change from 0.5 to 4 units in the violet-to-red peak intensity ratio (V/R), pointing out a signal of asymmetries in the disk. This was the first attempt to study the disk activity in π Aquarii, even though the stellar and disk parameters were undetermined.

Several authors (Nordh & Olofsson 1974; Underhill et al. 1979; Snow 1981; Goraya 1985; Theodossiou 1985; Fabregat & Reglero 1990; Zorec & Briot 1991) have studied this star over the years, trying to establish the fundamental parameters by using different observational techniques, such as multi-color photometry, interferometry, polarimetry, and high-resolution spectroscopy. In summary, based on all of these works, the effective temperature, T_eff, ranges between 22 500 and 30 000 K; the effective gravity, log g, is between 3.3 and 4.0 dex; the terminal velocity of the wind, v_∞, has been calculated as ~1450 km s⁻¹; and the mass-loss rate is $\dot{\text{M}}$ ~ 2.61 × 10⁻⁹ M_⊙ yr⁻¹ (Snow 1981).

McDavid (1986) used polarimetric data and found that the previous active phase started approximately in 1950, reaching its maximum around 1985, where the polarization level was higher than 1%. After that, as the emission line decreased, so did the brightness of the central star, reaching a minimum value in 1996 and extending its quasi-normal phase until 2002. Hanuschik et al. (1996) proposed that π Aquarii has a companion due to the broad and complex Hα line profile observed. Bjorkman et al. (2002) performed long-term observations of the star, obtaining information about the behavior of the disk phases of π Aquarii. They calculated new fundamental parameters of the binary system for the first time by fitting the photometric data with Kurucz local temperature equilibrium (LTE) stellar atmosphere models (Kurucz 1994). The best fit they obtained had the following values: M₁ = (11 ± 1.5) M_⊙, T_eff = (24 000 ± 1000) K, and A_V = (0.15 ± 0.03) mag. They used a distance of ${340}_{-70}^{+105}$ pc (ESA 1997), log (L_bol/L_⊙) = (4.1 ± 0.3) dex, and a R_⋆ = (6.1 ± 2.5) R_⊙. They also suggested an inclination angle in the range 50° ≤ i ≤ 75° based on previous studies on the FeII 5317 Å double-peaked emission line and taking into account that the central depression observed in Balmer lines between 1989 and 1995 is not as deep as it appears if seen edge-on. By examining the periodical radial velocity (RV) variation of the Hα emission line, they also concluded that if π Aquarii is a spectroscopic binary, a gas envelope should surround the companion. Considering this fact, they derived the following binary parameters: M₁ sin³ i = 12.4 M_⊙, M₂ sin i³ = 2.0 M_⊙, a = 0.96 arcsin i AU, and M₂ = 0.16 M₁.

Zharikov et al. (2013) recollected data for 40 orbital cycles from 2004 to 2013 and analyzed variations of double-peak asymmetries in Hα profiles. They searched for a periodicity in the V/R, finding 84.1 days, which agrees with the orbital period of the binary companion of π Aquarii (Bjorkman et al. 2002). In addition, they used evolutionary models to calculate the stellar parameters of the binary system. They obtained M₁ = (14 ± 1) M_⊙, M₂ = (2.31 ± 16) M_⊙, T_eff = (24 000 ± 1000) K, log L/L_⊙ = (4.7 ± 1) dex, R₁ = (13 ± 1.4) R_⊙, mass ratio q ≥ 0.165, 65° ≤ i ≤ 85°, and a = 0.96 AU, with a distance of 740 ± 90 pc, which is more than twice the previously reported value. They studied the disk’s structure by employing Doppler tomography, which showed a non-uniform structure. The structure they observed consists of a ring with inner and outer radii (in velocity units) of 450 and 200 km s⁻¹ and an extended stable region. The brightness of this extended region presented an S-wave that appeared similar to a one-armed spiral density pattern structure in phase with the secondary. In addition, the computation of the disk radius of ~65 R_⊙ (~0.33 AU) was estimated, assuming the Keplerian motion of a particle in a circular orbit at the outer disk.

Nazé et al. (2019b) collected data between 2013 and 2019 from the Be Star Spectra (BeSS) database¹ (Neiner et al. 2011) and found that the emission line profile seen in Hα almost disappears in January 2014. Then, the emission resembles the active phase seen from 1950 to 1990. A striking result was that the extended stable region (observed by Zharikov et al.) is no longer present after 2016, proving a change in the density distribution of the disk. In addition, the disk’s inner edge was closer (in velocity) to the central star, suggesting an increased disk size.

Observations from space-based photometric data have shown that pulsation is a common feature of classical Be stars, with these stars primarily pulsating in low-order gravity g modes, where gravity serves as the restoring force. Internal gravity waves transport angular momentum from the stellar interior to the envelope, causing an increment in the stellar rotation. This could create the conditions necessary to initiate an outburst, during which mass and angular momentum are transferred to the disk (Labadie-Bartz et al. 2022). Nazé et al. (2020) analyzes the photometric data of π Aquarii, revealing low-amplitude, coherent variability at high frequencies, suggesting p modes (pressure acts as the restoring force). This is atypical for hot Be stars because their internal structure and high rotation rates create conditions that do not favor the formation of stable pressure gradients in the outer layers of the star. However, other classical Be stars are also observed to have p mode pulsations (Walker et al. 2005; Huat et al. 2009; Labadie-Bartz et al. 2020).

Using observations with the spatial telescope XMM-Newton, Nazé et al. (2017) studied the variable thermal X-ray emission of π Aquarii, classifying it as a γ Cas type. These are a subgroup of classical Be stars with hard (more than 5–10 keV) and moderately strong continuous thermal X-ray flux. It has been suggested that the X-ray emission originates from accretion onto a secondary star. To determine the source of the unusual X-ray emission, simultaneous optical and X-ray observations were performed by Nazé et al. (2019a) from April to December 2018. In their comparison, the equivalent width (EW) of the Hα line increased from −22 to −24.5 Å, in contrast to the X-ray fluxes, which do not display correlation with time nor the EW of the Hα line. They note that although Hα has undergone dramatic changes, no significant changes were observed in the X-ray emission level, which is contrary to what has been found in other γ Cas stars. Moreover, they concluded that the behavior of π Aquarii points to a missing ingredient in the scenarios they considered.

Regarding the companion of π Aquarii, Tsujimoto et al. (2023) proposed that the secondary star is a non-magnetic or magnetic accreting white dwarf (WD). This information was obtained through two methods. (i) By comparing effective temperature and stellar luminosity to evolutionary models, they obtained an evolutionary mass of 10 M_⊙ (for the Be star). (ii) By giving a range of possible mass for the Be star (9–15 M_⊙) and for the inclination angle (50 to 75°), they obtained a dynamical mass for the secondary between ~0.8 and 1.4 M_⊙. Other works support the origin of X-rays in γ Cas stars due to accretion onto WD companions (Gies et al. 2023; Klement et al. 2024).

In this work, we study the line profile variability of several spectral lines identified in the spectra of π Aquarii in order to characterize the behavior of the star and the disk. We measure the EW, V/R, and double-peak separation (DPS) of the emission lines. For this study, we used spectroscopy data from the Be Stars Observation Survey (BeSOS) and BeSS databases (described in the next section). We derived a period from the V/R. We also obtained the stellar parameters under consideration of the oblateness of the star (limb and gravity darkening) and characterized the changes in the disk (density, inclination, etc.) by using the BeATLAS grid and Bayesian methods.

2 Observational data

The BeSOS database² (Arcos et al. 2018) contains optical spectra of southern Be stars, observed with an echelle spectrograph of R ~ 17 000. The catalog is complete up to magnitude V ~ 11. For each star, the catalog gathers information about its parameters, such as coordinates, magnitude, temperature, gravity, projected rotation velocity, and spectral type, among others. A total of 7 spectra of π Aquarii are available from the BeSOS database, where two were observed on November 14, 2012; three on July 24, 2013; and the last two on October 24, 2015. The spectra were averaged each night, obtaining three spectra with an effective wavelength range of 4200 to 7200 Å.

The Be Star Spectra database (Neiner et al. 2011) provides a comprehensive catalog of classical Be stars, Herbig Ae/Be stars, and B[e] supergiants. This database compiles spectra collected by both professional and amateur astronomers. Specifically, for π Aquarii, we have a total of 478 spectra from 2001 to 2024, covering only the Hα wavelength region and a total of 16 spectra covering the optical region. These include a spectrum from 2001 and at least one spectrum per year from 2008 to 2024. The wavelength range varies for each spectrum, as detailed in Table A.1.

Figure 1 contains all Hα emission line profiles plotted to appreciate the variation over the years, where epochs are sorted from the top to the bottom. This figure shows the high variation of this star over the years. Lighter colors (yellow) indicate higher intensity values. A great difference is observed between 2015 and 2024, where the intensity of the line increases considerably compared to the previous years. Details about this observed variation are presented in the results section under the EW, V/R, and DPS variability.

Fig. 1 Dynamical plot of Hα emission line profiles observed between 2001 and 2024.

3 Method

Several spectral lines were identified in the available spectra for π Aquarii. Some lines constantly appear in absorption or emission, while others exhibit both behaviors. Table 1 provides a detailed overview of these identified spectral lines, including descriptions of weak and strong lines, as well as variations such as double-peaked profiles. The following sections outline the methodology used to determine the stellar parameters and measurements on several lines and derive the disk parameters.

3.1 Determination of stellar parameters

The spectra of Be stars in the active phase are contaminated by emissions from the ionized disk, affecting not only hydrogen lines but also helium and other metallic lines. Therefore, to determine the stellar parameters, we chose the deepest absorption line, HeI 4471 Å (July 23, 2013), among the optical data available covering that range. The HeI 4471 Å is one of the most intense line profiles within all helium lines in the optical spectrum of O and B stars and is mainly constituted by a blend of a triplet 2³P–4³D and a forbidden 2³P–4³F transitions. As the radiation field of OB stars is not intense enough to populate the upper level of this transition (Osterbrock 1989), this line is, in principle, less perturbed by the contribution of the circumstellar environment. We used a new Python version of the ZPEKTR code (Levenhagen et al. 2011; Levenhagen 2014; Levenhagen et al. 2024), which adopts a grid of non-LTE spectra as input, built up with the SYNSPEC code (Hubeny & Lanz 1995), and based on the stellar atmospheres models calculated with the TLUSTY code (Hubeny 1988). The ZPEKTR code receives as input the radius R/R⊙, the mass M/M_⊙ and the effective temperature T_eff of the unperturbed, non-rotating spherical star and also the stellar rotation rate v/v_c and its inclination angle i. The star is divided into hundreds of thousands of grid elements, which are rotated at the specified input velocity rate and geometrically deformed and tilted. Each mesh element on the stellar surface has its particular radius value, computed with a Newton-Raphson iteration scheme, and its particular set of local physical parameters (T_eff, log g, v), which are evaluated obeying the classical von Zeipel’s relation with index β = 0.25. In the next step, the code performs the numerical inter-polation, within the mesh, of all local specific intensities using the non-LTE grid of TLUSTY/SYNSPEC models, accounting for a desired limb-darkening relation, and Doppler shifts these intensities accordingly. The output averaged flux spectrum is obtained by the numerical integration of all specific intensities emerging from the visible surface elements and brings the signatures of the effects of gravity darkening, limb darkening, and oblateness of the star. For this work, we used the linear Limb-darkening law (Klinglesmith & Sobieski 1970) and the Espinosa-Lara prescription for gravity darkening (Espinosa Lara & Rieutord 2011; Espinosa Lara 2014). We created a grid of 11 294 models with ZPEKTR code to seek the best representation of the stellar parameters. The range of parameters explored is described in Table 2, where the values are distinguished between the pole (subscript “p”) and the equator (subscript “eq”).

To select the best models, we conducted a Chi-square statistical test comparing the synthetic line generated using the ZPEKTR code with the observed spectral line. This analysis focused on the spectral range encompassing the HeI 4471 Å line. The determination of the “best fit” models was based on examining the change in the slope of the Chi-square values plotted against the number of models considered.

3.2 Estimation of equivalent widths

Before measuring the EW, all spectra were normalized to the continuum. We used different methods depending on whether the observed line is in emission or absorption. In the case of emission lines, that is, Balmer and some HeI lines (λ 5015, 5047, 5875, and 6678 Å) and some FeII lines (λ 4232, 4584, 5169, 5198, 5235, 5276, 5316, and 5363 Å), we subtracted the contribution of the stellar photosphere by using the result obtained with ZPEKTR. After removing the stellar photosphere contribution, the EW was computed considering the zero moments (M₀). This method quantifies the width (x_i) of a hypothetical rectangular line (of flux F_i) with a unit amplitude that has the same integrated area as the observed spectral line³ (Balona et al. 1996). We note that while subtracting the stellar contribution to the emission line, the continuum is set to zero, and then we do not subtract the unit from the resulting value. Despite this, we added a negative sign to follow the formal definition of EW for emission lines; that is, the more negative the EW, the stronger the line emission. To estimate the uncertainties of the EWs, the Bootstrap re-sampling method (Efron 1979) was implemented. The method constructs hypothetical data sets derived from the observed data, picking random points with replacement. Any data points may occur one or more times or may not be selected. The data points are size N, which is equivalent to the size of the observed data. The process was repeated 1000 times. The final value is given in the form: EW ± σ*, where σ* is the standard deviation of the entire bootstrapped data set.

3.3 Violet to red peak intensity ratio and double-peak separation

For all emission lines shown in Table 1, we calculated the V/R, which is the ratio between the violet and red emission peaks of circumstellar emission lines (Štefl et al. 2009; Carciofi et al. 2009). The DPS is the distance between the wavelength points of the violet and red emission peaks. The Hα line profiles present a complex geometry over time. Some profiles present a flat shape in the peak intensity, while others present a third peak in emission. To achieve this calculation, we classified the observed emission lines in three different shapes (cases) described as follows:

• Case 1: This case is for well-distinguished DPS profiles. An example of this type of Hα profile is shown in the left panel of Fig. 2. The method was to fit two Gaussian profiles, one for each peak. The center of each Gaussian gives the violet and red peak values, respectively.

• Case 2: Unlike Case 1, in this profile, the central absorption is not seen clearly, and it looks more similar to a third peak but with the violet and red peaks still distinguished on each side. The method used here was to locate the geometrical center of the line and the maximum emission at this center. The code then looks for the maximum intensities from each side of the center. The middle panel in Fig. 2 shows an example of this method.

• Case 3: This case is similar to Case 2 but has a flat emission line, making it difficult for the code to distinguish between the two peaks. The intensity peaks were selected by visual inspection between the extremes from each side of the geometrical center of the line (see the right panel in Fig. 2). We note that combining automatically computed values with those obtained by visual inspection can lead to inconsistencies and systematic biases. In particular, this could cause a larger scatter in our results. To diminish this possibility, we emphasized that higher peaks on each side were to be selected, and in cases where all peaks (when there are more than two) have similar intensities, the extremes of the peak emission were selected.

Fig. 2 Method used to obtain V/R and DPS for different shape profiles (solid black line). Case 1 (left panel): distinguished DPS profile. The Gaussian fits are represented in pink and red for the violet and red peaks, respectively. The dashed lines indicate the Gaussian fit’s center and the values for both peak intensities. Values for this example are in the legend of the plot. Case 2 (middle panel): the code selects the maximum intensities from each side of the center (dashed vertical lines). Case 3 (right panel): flat emission line. For these cases, the peaks were selected by visual inspection from each side of the geometrical center of the line (dashed vertical lines).

3.4 Errors for the violet to red peak intensity ratio and double-peak separation

The spectra used in this work come from different instruments with different spectral resolutions. Table A.1 contain a list of all instruments and spectral resolutions for the available data. To account for the fitting errors on the measurements of V/R and DPS considering the different spectral resolutions, we used an open module from Python⁴ to construct a wavelength grid by specifying the resolution (R = 20 000) and then rebinning the spectrum computing the mean within each bin, to the lower resolution in our sample (R=5800). Then, we measured V/R and DPS in both profiles and calculated the percentage error in Case 1 (and 2) for V/R and DPS, errors are 3% (and 2%) and 12% (and 4%), respectively. For case 3, we are considering the maximum error estimation found between both cases. In summary, the uncertainties have accuracies of 3% and 12% for V/R and DPS, respectively.

3.5 Periods

The high-cadence Hα line emission data for π Aquarii presents a significant opportunity to analyze the periodicity of the V/R emission line variation. By employing the astropy library’s implementation of the Lomb–Scargle periodogram (Lomb 1976; Scargle 1982), a sophisticated algorithm designed for the detection and characterization of periodic signals, we have successfully determined the period of these V/R shift variations. A key advantage of the Lomb-Scargle periodogram is its capacity to significantly reduce frequency aliasing, a phenomenon where higher frequency signals are misrepresented as lower frequencies due to inadequate sampling rates. This capability is crucial for accurately identifying periodic components, particularly those with frequencies exceeding the Nyquist frequency (Nyquist 1928). In addition to the Lomb-Scargle periodogram, we use the date-compensated discrete Fourier transform (Ferraz-Mello 1981) that offers a robust method for analyzing unevenly spaced data. The date-compensated discrete Fourier transform is specifically adapted to handle irregular temporal sampling by employing the concept of function space projection to perform a Fourier transform.

3.6 Determination of disk parameters

To estimate the disk parameters of π Aquarii, we used a grid of models created using the 3D Monte Carlo radiative transfer code HDUST (Carciofi & Bjorkman 2006, 2008) to model Hα lines and the spectral energy distribution (SED). This code solves the radiative and statistical equilibrium equations for hydrogen, considering non-LTE and different types of opacity and considering fast rotation for the central star; therefore, it assumes gravity darkening and oblateness effects. To model the circumstellar gaseous disk was used the steady-state power-law approximation of a viscous decretion disk (VDD) model (Lee et al. 1991; Bjorkman 1997): the density distribution falls off as a power law radially (r) and a Gaussian distribution in hydrostatic equilibrium axisymmetric respect to the mid-plane of the disk in the vertical direction (z). In this approximation, the volumetric density distribution can be written as $$\rho (r,z)={\rho }_0{\left(\frac{r}{R_{eq}}\right)}^{-n}\exp \left(\frac{-z^2}{2H^2}\right),$$(1)

where ρ₀ is the disk base density, R_eq is the equatorial radius of the central star, n is the volume density distribution index and H is the height scale that sets the flaring of Be stars and is given by $$H=H_0{\left(\frac{r}{R_{eq}}\right)}^{\beta }.$$(2)

For an isothermal disk, β=3/2 and the constant H₀ is defined as $$H_0={\left(\frac{2R_{eq}^{3}k{T}_0}{GM{\mu }_0{m}_H}\right)}^{1/2},$$(3)

where M is the stellar mass, G is the gravitational constant, m_H is the mass of a hydrogen atom, k is the Boltzmann constant, μ₀ is the mean molecular weight of the gas and T₀ is an isothermal temperature used only to fix the initial vertical structure of the disk.

With a vertically integrated disk density ρ(r, z), we can obtain the surface mass density Σ(r) as follows: $$\mathrm{\Sigma }(r)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\rho (r,z)dz.$$(4)

The relation between the volume and surface mass densities at the base of the disk Σ₀ is given by $${\rho }_0={\mathrm{\Sigma }}_0\sqrt{\frac{GM}{2\pi c_s^2{R}_{eq}^3}}.$$(5)

As input the code takes three disk parameters: the number density⁵ n₀, the radial exponent n, and the disk truncation radius R_d. The isothermal temperature was fixed at T₀ = 0.72 T_eff (Correia Mota 2019), and the mean molecular weight of hydrogen was μ = 0.6.

A source for the central star needed to be defined. We selected a parametric rigid rotator providing the following parameters: M, R_p, R_eq/R_p, L, log g, V_rot, and T_eff. The stellar atmosphere models of Kurucz (1994) are used for the stellar spectrum and solar metallicity (Z = 0.014) was used. The limb darkening models from Claret (2000) are adopted to outline the photospheric specific intensity as a function of direction.

The HDUST code was used to create a grid of models called BeAtlas (Rubio et al. 2023). Our version of this grid contains a total of 13 800 models, which include Balmer lines and SED in the range between 0.10 μm and 2.45 μm. The age adopted for the grid models is fixed and has been parameterized using the fraction of hydrogen in the core X_c = 0.3, which corresponds to the end of the main sequence, where most Be stars are concentrated.

The model grid adopts the normalized value (between 0 and 1) for the base surface density of the disk Σ₀ during modeling for two reasons: the lower and upper limits are separated by a factor of ~10², which makes it difficult to explore the parameter space and the fact that more massive stars have much denser discs when compared to less massive stars (Rímulo et al. 2018).

The parameters contained in the model grid are 11 values of masses (M), which correspond to spectral types between B0.5 and B8 (Townsend et al. 2004), 5 values of oblateness (R_eq/R_p) which corresponds to different rotation rates of the star, the normalized base surface density of the disk (Σ₀), 4 values of volume density distribution index (n) which includes a realistic scenario of the disk, with values that indicate a dissipation (n = 3), stationary (n = 3.5) and growth (n = 4.0 and n = 4.5) of the disk (Haubois et al. 2012), and 10 values of the cosine of the inclination angle of the disk concerning the observer’s line of sight, cos(i), the value equal to 1 corresponds a pole-on view while the value equal to 0 an edge-on view. A range of each parameter is shown in Table 3.

In addition, through the modeled parameters, the use of Geneva stellar evolution models (Georgy et al. 2013; Granada et al. 2013) we can estimate other stellar parameters such as luminosity (log L) and polar radius (R_p). The effective temperature at the pole (T_p) using Eq. (4.21) from Cranmer (1996), polar surface gravity (log g_p) from its relationship with the mass and polar radius, the fraction of critical angular velocity (ω) using Eq. (9) from Ekström et al. (2008) and gravity darkening exponent (β_GD) through the model presented by Espinosa Lara & Rieutord (2011), where the exponent depends on the rotation rate W.

To determine the model that best fits our observations, we used the Bayesian code EMCEE⁶ (Foreman-Mackey et al. 2019), which uses a variation of the Markov chain Monte Carlo (MCMC) algorithm from Goodman & Weare (2010). This code is an ensemble sampler. It uses the concept of “walkers” to explore the multidimensional parameter space. Numerous of these walkers, defined by the user, are started together and, throughout the iterations, start exploring the parameter space independently (but connected), creating large chains that tend to converge on one value or more (in cases where a bi-modal is observed). As a result, we obtain the probability density functions (PDFs) for each parameter of our grid, which we use to estimate the errors and observe possible correlations between them. Among the advantages of using EMCEE code, we can mention i) high efficiency when modeling multidimensional parameter space and ii) parallelization of operations, reducing computation time.

According to Bayesian formalism, we must define the likelihood function and a priori distribution. We used the likelihood function, ln P(D|θ, I), of the χ²-distribution type given by $$lnP(D|\theta ,I)\propto -\frac{1}{2}{\chi }^2,$$(6)

where χ² is defined as $${\chi }^2=\sum _{i=1}^N\frac{{(F_{obs,i}-F_{mod,i})}^2}{{\sigma }_{obs,i}{}^2}.$$(7)

The logarithm of the prior distribution, ln P(θ, I), corresponds to the knowledge available about each of the parameters in our grid or any other relevant information that helps the modeling process. To follow the same range of parameters used in the stellar parameters determination, we have adopted a non-informative prior distribution for the mass range between 10 and 14.6 solar masses and inclination of the disk between 50 and 75 degrees, which correspond in terms of cos(i) to a range between 0.25 and 0.65 according to the information available in the literature (see the introduction). For the other 3 parameters, we also adopt a uniform prior distribution, but that covers all the values available in our grid. Another parameter we adopted as a prior distribution was the v sin i; in this case, we use a Gaussian distribution (similar to Eq. (7)) with mean 227 km s⁻¹ and standard deviation 78 km s⁻¹ from Solar et al. (2022). For SED modeling, we include two other parameters: E(B-V) and the parallax, both using a Gaussian distribution according to the values obtained by (Fitzpatrick & Massa 1999) and Gaia DR3 (Gaia Collaboration 2023), respectively.

We modeled a total of five Hα lines obtained at different years (i.e., 2001, 2011, 2014, 2018, and 2022) and SED using photometry data available in the Virtual Observatory VOSA (Bayo et al. 2008) from missions and observers in ultraviolet, including the IUE mission from Benvenuti (1983), visible; UBV photometry from Mermilliod (2006); Stroemgren V filter from Paunzen (2015); Tycho-2 Catalogue from Høg et al. (2000); Hipparcos catalog (ESA 1997); and Gaia DR3 (Gaia Collaboration 2016), and in the near-infrared the 2MASS catalog from Skrutskie et al. (2006). The selected range of 0.14–2.4 μm, corresponds to the spectral coverage of our grid.

A total of 150 walkers were defined for all the modeling. To check the convergence of the walkers, we followed the recommendations given in the EMCEE documentation⁷, where the number of iterations reaches the value corresponding to 50 times the correlation time and shows a variation of less than 1% in subsequent iterations.

Fig. 3 Observed spectrum of the photospheric line HeI 4471 Å (in black). The red line shows the best fit of the ZPEKTR models, and the gray band represents the flux values of the 67 best models.

4 Results

4.1 Stellar parameters obtained

We created 11 294 models with ZPEKTR, where the best models are 67. We obtained the average and standard deviation from the output parameters. The best-fit models are shown in Fig. 3, where the gray band indicates the range of the best 67 models. The stellar parameters derived from these models are presented in Table 4. The uncertainties were calculated as the standard deviations from the 67 models compared to the best-fit model.

Since W is not provided by ZPEKTR, it can be calculated using Eq. (11) from Rivinius et al. (2013). Applying standard error propagation theory, we obtained a final value of W = 0.55 ± 0.01.

4.2 Obtained disk parameters

The observed profiles of the Hα lines are normalized in flux. For this reason, we did not include parallax and E(B − V) as additional parameters in the modeling since, for an isolated line, we can disregard any significant changes in its profile. The profiles were interpolated to the same resolution as the models contained in the grid. In addition, we used the velocity space, centered on the Hα line, with a range between −1100 and +1100 km s⁻¹, sufficient to show the line profile and its wings.

Figure B.1 (upper right inset) shows the results obtained for the spectrum observed in 2001. We can note that the PDFs are defined for mass, oblateness, cos i, and normalized Σ₀, the latter with a smaller uncertainty than the others. The n index showed an irregular distribution, almost flat. The line profile observed was well-fitted, both in the central region and on the wings. Using the Spearman coefficients⁸, we can also obtain information about the correlation (positive values) or anticorrelation (negative values) between each pair of parameters. The higher the value and the intensity of the color, the stronger the correlation between the parameters. We can see strong correlations between the normalized base density Σ₀ and the n index.

Figure B.2 shows the results obtained for the spectrum observed in 2011. On this epoch, we can observe an intense double-peak profile with V/R ~ 1. The fit obtained by the model was well-defined, as were the PDFs for each parameter. Similar to the results of 2001, an intense correlation between the disk parameters can also be observed.

Figure B.3 shows the results obtained for the spectrum observed in 2014, where we can note a bimodality in the PDFs of three parameters, with the presence of very pronounced peaks, except for Σ₀ and cos i, this can indicate that both solutions are probably valid. The correlation/anticorrelation between the parameters is complex, probably due to bimodality. The modeled line profile reasonably agrees with the observed profile, especially on the left wing, where contamination appears to increase the flux in this region.

Figures B.4 and B.5 show the results for the spectra observed in 2018 and 2022, respectively. In both cases, the observed line profiles were not well-fitted. This is due to the complexity of these lines, which have a width greater than 500 km s⁻¹ and asymmetrical V/R values. Our version of the grid has limitations, including the lack of models with asymmetric line profiles caused by inhomogeneities in the disk and the non-inclusion of physical processes, such as electronic scattering, which would broaden the line profiles. Although the PDFs are well-defined, with pronounced peaks, the values obtained for the parameters are not reliable.

Figure B.6 shows the results of the SED modeling. Although the photometric data were obtained at different epochs and could be affected by the variability caused by the disk formation/dissipation process, we obtained a good fit. One possible explanation may be related to the coverage of the grid between the ultraviolet and the near-infrared since, in this region, the star’s parameters are better defined than the parameters of the disk. Therefore, they are less affected by photometric variations coming from the disk. This is evident from the poorly defined PDFs of the disk parameters. Both E(B-V) and parallax, included in the SED modeling, showed well-defined PDFs. Tables 5 and 6 contain stellar and disk parameters, respectively, obtained as a result of the modeling of the Hα line observed on different dates and the SED.

4.3 Measurements obtained over Ha emission line: EW, V/R and DPS

4.3.1 Equivalent width variability

Figure 4 shows the evolution of Hα EWs, V/R, and DPS from 2001 to 2024. The EW values of Hα were obtained for all observation epochs and corrected for photospheric absorption using synthetic spectra from ZPEKTR. These values are available in an online spreadsheet. Thus, the Hα EWs values represent ‘pure’ emission, and therefore, the growth of the EW curve will indicate the disk growth, as well as the decrease of the EW curve, which will indicate a disk dissipation phase. We note to the reader that all results in this work do not consider the effects of the continuum in the calculated EW. This is because most spectra are already normalized to the continuum in the BeSS and BeSOS databases. Possible effects are described in the discussion section. The sharp increase in recent years may be related to episodes of enhanced mass-loss rates or outflows. Gray vertical lines in Fig. 4 indicate five significant epochs where noticeable changes in EW occur, providing reference points to help explain the time variability: MJD⁹ 52265 (December 22, 2001), 55837 (October 03, 2011), 56959 (October 29, 2014), 58400 (October 09, 2018), and 59908 (November 25, 2022).

Due to the lack of earlier data in this study, no conclusions can be drawn about the growth or dissipation of the disk before 2001. However, Bjorkman et al. (2002) reported the lowest brightness level observed in π Aquarii in 1998, along with weak Hα emission, which they described as a ‘quasi-normal star’ phase between 1996 and 2001. Between December 2001 and October 2011, the EW grew slowly and fluctuated around −10 Å. Then, from October 2011, the EW decreased, reaching a minimum documented value close to −2 Å and making an inflection point (in October 2014 and marked with a red dashed vertical line) in the behavior curve where it begins to increase its value steadily over time until it reaches a value of −28.7 Å in October 2018 (at the fourth gray vertical line). Subsequently, the EW fluctuates around this value until it increases drastically, reaching a value of −36.1Å in November 2022 (at the fifth gray vertical line), which is more than 18 times the value at the inflection point. This suggests that the growth of the circumstellar envelope has continued without interruption since October 2014; an outflow or enhanced mass-loss rate occurred in recent years (November 2022). After this, the EW started to decrease rapidly.

Table 7 displays the EW ratio between the Balmer lines and Hα for the five selected observation dates. The results show that the ratio amplitude increases toward the Hγ line, with similar EW behavior observed for all the Balmer lines.

Fig. 4 Equivalent width, V/R, and DPS of Hα plotted against time in MJD. Gray vertical lines denote five specific dates with significant changes in EW, which are discussed in the text. A horizontal line at V/R = 1 is included to highlight the asymmetry of Hα.

4.3.2 Violet to red peak intensity ratio variability

The V/R presents periodical fluctuations with intensity values between 0.69 and 1.24 over time. Despite this, these variations have some differences depending on the epochs. For example, between October 2011 and 2014 (second and third gray vertical lines), in the zones with a high concentration of observations, the values describe an M-type shape (typical of a cyclical variation) around one. We note that before the EW plot’s inflection point (October 2014), the V/R intensity ratio reached a higher and deeper value than the previous months. After the inflection point, the lowest values are less deep, and the highest values are more intense, becoming asymmetric and biased toward the violet peak. Since October 2014 (third gray vertical line¹⁰), the V/R ratio has become confusing as the disk grows. This can be a result of the precession due to the one-armed density wave presented in the disk. These values could also result from losing the symmetrical shape of the Hα emission line; as can be seen in the DPS curve, values are lower, indicating a more extended region of the emitting area.

The Lomb-Scargle periodogram shown in Fig. 5 has been computed for the V/R of Hα. The left panel shows the power spectrum. The maximum amplitude corresponds to 81.75 days for the date-compensated discrete Fourier transform periodogram and 81.94 days for the Lomb-Scargle periodogram. The other maximum frequency corresponds to 69 days. The right panel displays the phase of the V/R from 0 to 1, calculated with the best period of 84.75 days, which exhibits a sinusoidal shape. Zharikov et al. (2013) searched the period of the V/R using a Discrete Fourier Transform method, and the dominant frequency they found in the power spectrum was 0.01187 cycles/days, corresponding to a period of 84.2 days, which agrees with our results. Our result in this work agrees with the period of the orbital system of 84.1 days found by Bjorkman et al. (2002).

4.3.3 Double-peak separation variability

By measuring the DPS, information about the rotational velocity, size of the emitting region, inclination, temperature and density distribution, and disk dynamics can be inferred. Because DPS is primarily a result of Doppler shifts due to the Keplerian rotation of the disk, faster rotation leads to a larger separation, implying that the emission originates closer to the star in the inner disk (Huang 1972). We note that this is also sensitive to the inclination angle of the disk relative to the observer, whereas, for edge-on disks, the separation is wider due to the higher radial velocity components visible to the observer. For π Aqr, the DPS curve goes from a high value, close to 530 km/s, to a lower value of ~100 km/s, in all analyzed periods. In detail, from December 2001, the DPS continuously drops its value until October 2011 (second gray vertical line), where it remains almost constant at around 300 km/s. After the inflection point, the DPS decreases continuously until it reaches an almost constant value around ~150 km/s. Therefore, based on Hα measurements, the emitting region increased ~18 times its value since the inflection point and decreased by nearly twice its rotational velocity.

Fig. 5 Top panel: power spectrum for the V/R presented in units of periods. The date-compensated discrete Fourier transform periodogram is shown as a black solid line, while the Lomb–Scargle periodogram is represented as a blue solid line. Both periodograms have been normalized by their respective means to enhance visibility. Bottom panel: phase of the V/R computed with the period of 84.75 days.

4.4 Other lines

As we mentioned in the observational data section, not all observations downloaded from BeSS and BeSOS contain the full optical range. Therefore, it is not possible to compare the other studied lines from Table 1 with the exact five selected epochs of Hα. Hence, from the 16 spectra covering the optical region, we chose 14 dates to compare simultaneously the shape of the lines (the remaining two were not chosen because they were at consecutive dates where no considerable variation was seen). The selected elements presented emission at least once: helium, silicon, and iron, as well as Hβ and Hγ. From Hα line, we infer that all chosen dates are in a disk phase. The line profiles are presented in velocity, and the intensity was normalized to the continuum. We remember that the lines are all corrected by the contribution of the atmosphere. The related figures are displayed in the Appendix C.1 for Balmer, Appendix C.2 for helium, Appendix C.3 for silicon, and Appendix C.4 for iron. In the following, we compare the shape of the line profiles in detail concerning the Balmer lines. We note that three of the five disk phases discussed above coincide with the observation dates chosen here for the other lines (October 09, 2018, and November 25, 2022).

4.4.1 Hβ and Hγ hydrogen lines

Figure C.1 shows Hα, Hβ, and Hγ line profiles for different disk phases. The vertical axis is different for Hα due to its strong emission. From this figure, we see that when the emission is presented, a DPS is always observed; however, the V/R intensity ratio is not necessarily equal for all Balmer lines, and the intensity of the emission is not always descending with respect to higher orders in the series (higher-order Balmer lines require higher densities and temperatures to emit significantly; thus they primarily trace the inner, denser disk). For example, in October 2011 (which also coincides with one analyzed disk phase at the second gray vertical line), Hβ is less intense than Hγ. This repeats in November 2020 and September 2022. Moreover, on these three occasions, the DPS is also different, and it seems that Hβ does not follow the same symmetry as Hα and Hγ. This suggests that the regions of the disk contributing to the emission of these lines differ significantly in their optical depth and/or geometry. Following the Hα emission line shape and the results for the V/R intensity ratio, a small overdensity is presented in the disk for years until the inflection point occurs (October 2014), where the disk density increases and its distribution changes.

4.4.2 Helium lines

Figure C.2 shows three selected helium lines from the Table 1, which show emission lines features. These helium lines are HeI λ 5875, 6678, and 7065 Å. An empty box will appear if no data are available for the line on that date. In general emission lines show a shell profile with low intensity emission (≤ 1.1). We note after the inflection point, between 2016 and 2020, no emission (or a very low) is seen in λ 6678 and 7065 Å, but a persistent low emission is observed in 5875 Å. This line is more easily excited than the other two HeI lines, which are higher-energy transitions. Then, if the disk is dense and moderately heated but not extreme (we obtained a T_eff ~ 24 000 K, and if we consider that T_disk is between 60 and 72% of T_eff then T_disk is approximately between 14 000 and 17 000 K) the 5875 Å line, may appear in emission while the others do not. Denser and higher temperature disk regions are found close to the star; therefore, HeI emission lines are traced to the inner part of the disk or optically thick disk zones. The last observation date in Figure C.2, when Hα is higher in intensity and larger in EW, shows emission in the three HeI lines, with similar intensities and with a V/R < 1, as observed in Balmer lines.

4.4.3 Silicon and iron lines

Figure C.3 shows the silicon lines 6371 and 6347 Å and Fig. C.4 displays the iron lines 4232, 4584, 5018, 5169, 5198, 5235, 5276, 5316, 5363 and 6385 Å. The calculated values of the EW and DPS for spectral lines of both elements are available in the online excel file. Although the spectra of SiII lines are very noisy, a very low emission is distinguished between 2019 and 2022. The emission lines are asymmetric, showing a DPS. These low-ionisation lines form in cooler and denser regions of the disk compared to helium or hydrogen lines.

All FeII lines have very similar behavior except for the case of FeII 5018 Å, where an absorption line can be distinguished in the spectra between 2001 and 2013, with a small DPS profile in some observation epochs, similar to Hβ. Since emission appears in FeII 5316 Å, it is the most prominent compared to the other iron lines. Although the V/R is more similar to Hβ than helium lines, the shape of the lines shows a flat structure in the emission peak, with a third peak. In 2022, FeII λ 5018 and 6385 Å show opposite behavior in the V/R ratios compared to the other iron lines. Not all the FeII data exhibited V/R variations over the years. In some lines, no emission or an absorption scheme is observed. The most significant variation of the DPS is encountered for FeII 5018 Å, where the highest DPS is 585 km/s on 2010-10-26. Compared with Hα on the same day, it is almost twice as large. Following this, the DPS continuously drops its value to its minimum of 152 km/s. From FeII 5169 to 5363 Å, the DPS fluctuates from 412 to 205 km/s. The DPS of FeII 5198 Å continuously decreases from 372 to 204 km/s. All the other lines have variations between their maximum and minimum peaks. We remark that the highest value of DPS for each line occurred on 2015-10-31, with the highest DPS of 412 km/s for FeII 5276 Å. From FeII 5169 to 5276 Å the smallest DPS is observed in 2022-09-03. Any FeII line has a similar behavior decrease in the DPS compared to Hα, except for the FeII 5198 Å. The lowest DPS is found for FeII 6385 Å, equal to 87 km/s. This is the smallest DPS for all the lines in π Aquarii. For FeII 6369 and 6385 Å, the DPS is decaying similarly to Hα. In the discussion section below, we include Fig. 6 with the EW versus DPS of iron emission lines studied in this work.

Since the inflection point, almost all iron lines have been in emission. FeII emission is generally associated with intermediate-density regions of the disk, typically closer to the outer or mid-disk. Its transitions involve lower excitation energies and are less temperature dependent than SiII. By comparing the variations of the spectral lines of the other elements, the strong hydrogen and FeII emission lines, alongside weak SiII emission, suggest that the outer disk is dominant (the outer and cooler parts of the disk are contributing most of the observed line emission).

Fig. 6 Double-peak separation of iron lines compared to their EW for the penultimate observation date of Fig. C.4. This date corresponds to December 31, 2021. Red vertical dashed lines denote the DPS for the Balmer lines at the same date.

5 Discussion

5.1 Possible effects of not considering continuum contribution in the equivalent width determination

As we mentioned previously, most of the spectra used in this work are already normalized to the continuum in the BeSS and BeSOS databases. This is relevant because during an episode of mass ejection, the disk can grow, and the distribution of material in the disk can change. On the one hand, when the Be star is viewed edge-on, a thicker disk with more ejected material would cause a greater reduction in the continuum intensity because the material in the disk can block or scatter more of the star’s light. The continuum intensity may decrease, and the spectral features (e.g., emission or absorption lines) may become more prominent. On the other hand, when the system is viewed pole-on, the observer sees less of the disk material, and the scattering effects are reduced. The continuum observed would be closer to the star’s intrinsic continuum, and changes in the observed continuum due to mass ejection into the disk would be less pronounced. Emission lines may also be less affected because the amount of disk material in the observer’s direct line of sight is smaller. Therefore, the EW could be underestimated by the changes in the continuum, overall, if the Be star is seen from the edge-on. We note that we discuss and conclude relevant aspects based on the EW values obtained without taking the continuum into consideration. In particular, when we claim that the disk is dissipating (because the EW is lower), it may not be completely true since the effect described above.

5.2 Relation between helium, silicon, and iron lines concerning the disk phases

Many iron lines present emission since the inflection point of Hα. Several works in the literature are related to the study of iron lines (e.g., Hanuschik (1987); Ballereau et al. (1995); Arias et al. (2006)). The first work divides the shape of the line profile of iron lines into two classes: Class 1 are double-peaked and symmetrical emission lines, which describe the relationship between the width of the peaks concerning v sin i. Class 2 line profiles exhibit single asymmetric sharp peaks in a complex, broad emission structure. It also concludes that the emitting region of iron lines is inside Hα and Hβ formation region zone, close to the central star and that they are optically thin.

It is possible to calculate the emitting region for spectral lines using the Huang (1972) relation: $${\mathrm{R}}_{\mathrm{H}}={\left(\frac{2\mathrm{v}\sin \mathrm{i}}{\mathrm{D}\mathrm{P}\mathrm{S}}\right)}^2\mathrm{R},$$(8)

where v sin i is the equatorial projected rotational velocity of the central star, R is the stellar radius, and R_H is Huang’s radius for the emitting region. This relation was used to compute Huang’s radius for the Balmer lines on the penultimate observation date of Fig. C.4, which correspond to 2021-12-31, where v sin i = 271 km/s, taken from Table 4. We have for Hα: R_H = 7.0R, Hβ: R_H = 4.5R and Hγ: R_H = 2.7R. Looking at the FeII lines, the only ones outside the emitting region of Hα and Hβ are FeII 5018 and 6249 Å, with R_H equal to 12.7 and 7.9 R, respectively. We also note that FeII 5018 Å is the only FeII emission line showing a different shape profile (V/R) compared to the rest of the lines (see Fig. C.4). In Fig. 6, we show the variation between EW and DPS for 2021-12-31 of iron emission lines indicated with different colors and symbols. The plot also shows the DPS of Balmer lines to compare with. We chose this date because all FeII emission lines show a clear DPS and have spectra on that observation epoch. Looking at the figure, we can conclude that there is no correlation between the EWs of iron lines and the disk formation region. The DPS shows a wide variation range. Almost all iron lines are concentrated on the inner disk inside the Hβ emitting region.

Ballereau et al. (1995) asserted that iron lines must be optically thick, which can change the location of the formation region where FeII lines are formed. Arias et al. (2006) studied the optical depth regime of FeII λ 5198, 5235, 5265, 5276, 5317, and 5365 Å lines and the degree of expansion of their formation region simultaneously, finding out that these lines are optically thick. They are formed in circumstellar regions close to the central star. This is also confirmed by the fact that the source function of FeII lines rapidly decreases with radii, preventing the formation in the outer part of the disk. They calculated the extent of the emission region of FeII lines using the self-absorption-curve (SAC) method (Friedjung & Muratorio 1987), which allows the opacity effect to be explicit on the emitted radiation intensity. This is in contrast to Huang relations, which are only for optically thin lines and result in uncertainty in the outcome presented before.

5.3 Evolutionary status of π Aquarii

We study the evolutionary state of π Aquarii with the results obtained with the code ZPEKTR, using the evolutionary tracks from Georgy et al. (2013); Granada et al. (2013). This grid was made for a range of masses between 0.75 to 15 M_⊙, separated in 270 evolutionary tracks, for rotating stars at Z = 0.014, 0.006 and 0.002, the interpolated models were obtained from this link. Figure 7 shows the Hertzsprung-Russell diagram for stars with ω = 0.80 and Z = 0.014, where the solid lines indicate the evolution tracks with the corresponding initial mass. Zharikov et al. (2013) computed the mass using the stellar parameters; they calculated the M_V with system brightness V equals 4.85 mag, A_V = 0.15 mag, and HIPPARCOS distance of 340 pc, and using a T_eff = 24 000 K with a bolometric correction BC = −2.36 mag, they finally obtained log L = 4.02 L_⊙ that leads to an initial mass of 10.5 M_⊙. Using a shorter distance of 240 pc (van Leeuwen 2007), the initial mass changes to 9.5 M_⊙ with the models from Ekström et al. (2012).

Our result is consistent with an initial mass between 10 and 11 M_⊙, located within the main sequence. This is about halfway through the main sequence. This result agrees with the result of ZPEKTER. The dynamical mass found by Bjorkman et al. (2002) is between 12.5 and 17 M_⊙, calculated in a disk-less phase with a range of i from 65°−85°. Our stellar parameters were computed in 2013, when the EW of Hα dropped to a minimum, with a range of i between 50°−75°. With these new models, the initial mass of Zharikov et al. (2013) would change from 10 to 11 M_⊙. Tsujimoto et al. (2023) using T_eff = 24 000 ± 10 000 and a Gaia distance of 286 pc obtained a value of log L = 3.87 L_⊙, placing it in the evolutionary state with an initial mass close to 10.5 M_⊙, where the range they proposed between 9 and 15 M_⊙ is consistent, the reason for this wide range is to include the dynamical and evolutionary mass.

Several authors have calculated log g and v sin i of π Aquarii. Chauville et al. (2001) studying the HeI 4471 Å line determined T_eff and logg of 116 Be stars. For π Aquarii they found T_eff= 26 668 K and logg = 3.95 dex. Frémat et al. (2005), utilizing a code that considers the fast rotation of stars as rigid rotators and gravity darkening, found T_eff = 26 061 ± 736 K and log g = 3.92 ± 0.09 dex. Arcos et al. (2018) used photometry data fitted with Kurucz (Kurucz 1994) and TLUSTY (Hubeny & Lanz 1995) stellar atmosphere models, together with spectroscopy data, to find the parameters of a Be star sample. For our star the results were: T_eff = 24 500 ± 245 K, log g = 3.4 ± 0.03 dex and v sin i = 215 ± 4 km/s. Ahmed & Sigut (2017) using photospheric nitrogen abundances, computed T_eff = 26 061 K and log g = 3.7. We note that different authors’ temperatures are similar, but the effective gravities in all works are always lower than 4.

Rivinius et al. (2013) summarize the main properties of the rotation of Be stars. They derive an average value of $\overline{\mathrm{W}}=0.8$, based on spectroscopic (Frémat et al. 2005) and interferometric (Meilland et al. 2012) analyses. They note that this value does not depend on temperature or effective gravity; more importantly, the lowest W value for a B star to become a Be star is about 0.7. The W calculated in our work is lower than these values, about 22%. Looking at other works, Zorec et al. (2016) computed the true rotational velocities of 233 Be stars, fixed by the effects of gravitational dimming, binarity, and macroturbulence, and found that the mode of γ is equal to 0.65. The value of γ of π Aquarii in Table 4 agrees with this value. Cochetti et al. (2019) derived the mean values of γ = 0.75 and ω = 0.9 of 18 Be stars using interferometric observations. In comparison with the latter authors, it is found that the values of γ = 0.63 and ω = 0.82 calculated by ZPEKTR of π Aquarii are smaller.

Figure 8 shows the R_eq/R_p ratio versus projected velocity for different Be stars from van Belle (2012), ranging from III to IV in luminosity class, which is in agreement with the spectral classification of π Aquarii, which is classified as B1III-IVe by Slettebak (1982).

Fig. 7 Hertzsprung-Russell diagram for π Aquarii. The solid lines represent the evolutionary tracks of rotating stars with solar metallicity and ω = 0.8, as outlined by Georgy et al. (2013) and Granada et al. (2013), the dash magenta line has an initial M = 10.5 M⊙. The filled red circle denotes the findings from this study. The filled blue circle indicates the luminosity and temperature derived from Zharikov et al. (2013), calculated using varying distances. The filled gray circle corresponds to data from Tsujimoto et al. (2023). The filled green circle corresponds to data from Hohle et al. (2010).

5.4 Quantification of the disk changes under parametric model assumptions

We noticed that the profiles of this star range from shell type to DPS to triple-peaked to flat-topped profiles. This is typical behavior for a Be star in a binary system (Okazaki et al. 2002; Marr et al. 2021; Panoglou et al. 2016, 2018; Cyr et al. 2017). Panoglou et al. (2018) used a smoothed particle hydrodynamics (SPH) code to compute the configuration of Be discs in coplanar circular binary systems. The output from SPH consisted of a disk configuration of two spiral arms. This result was used as input into HDUST code to obtain the theoretical Hα emission lines view at different inclination angles. They proposed that triple-peaked profiles are an evolved state of flat-topped profiles. Moreover, they stated that the evolution with the orbital phase of the profile shape can be illustrated as a DPS followed by a flat-topped profile then a triple-peaked profile and again a flat-topped profile and then a DPS profile. Then, π Aquarii can be in such a coplanar or misaligned binary system (Moritani et al. 2013).

By assuming the VDD’s parametric model, we obtained good results in modeling the Hα emission line for the epochs 2001, 2011, and 2014. This last presented a bimodal distribution in all parameters except for the inclination angle. In 2001, the superficial density was lower and decayed faster, close to the star. In 2011 and 2014, the density decays slower with n exponent values n ~3.2. In the last two epochs, 2018 and 2022, the model fit the intensity of the profiles but not the shape. On one hand, in 2018 the emission line is irregular with a flat-topped peak, on the other hand, in 2022 the emission line presents a DPS profile with R≫ V. Also, the wings are not well reproduced, mainly because HDUST does not reproduce electron scattering (a process that widens the wings of the line). For these two epochs, we can see that the n exponent does not converge correctly, giving the lowest possible value of n within the grid (n = 3.0). Smaller values of n are needed to achieve a slower density decay near the disk to reproduce the observed line. The superficial density increases with the years, as can be observed visually as a growth in the shape of the emission line. The inclination angle ranges between ~60° and 77° (or 70° if we consider only the acceptable fits). The possible reason for such a difference in the change of the inclination angle is that the disk is precessing and therefore is in a misaligned binary system.

Cochetti et al. (2019) used HDUST to model infrared spectra obtained in June 2017 with FIRE spectrograph and obtained good fits with the following parameters. n = 3.5, ρ₀ = 10⁻¹¹ g cm⁻³ and i = 45°. Comparing these results with the parameters obtained by us at a similar epoch (2018-10-09) (see Table 6), we found no similarities for the angle of inclination and the value of n. These differences may be related to inaccuracies in the modeling of the Hα line, as mentioned in the previous paragraph, and also by the approximations for the mass, radius of the star, and rotation rate used to create the infrared models. As for the parameter ρ₀, calculated using Eq. (5), we found a value of 1.75 × 10⁻¹¹ g cm⁻³, slightly higher than that estimated by Cochetti et al. (2019).

Fig. 8 Oblateness versus projected velocity of Be stars, as presented by van Belle (2012). The solid circle represents π Aquarii, with data derived from Table 4. Different colors indicate the luminosity class of each star, with the red color representing the star Pi Aquarii.

5.5 Origin of X-ray emission in π Aquarii: Accretion onto the secondary or interaction in the Be disk

As mentioned in the introduction section, π Aquarii presents hard continuous thermal X-ray emission, which can be the product of accretion onto the secondary (most probably a magnetic or non-magnetic WD), wind or magnetic interaction between the star and its decretion disk (Smith et al. 2016; Rauw et al. 2022). Our five selected dates mark a change in the EW Hα variation; therefore, we interpret these changes as a disk phase. We can obtain an estimation of the emitting region by using Huang’s law (see Eq. (8)). From our measurements, the DPS of the Hα emission line varies significantly across the five epochs, ranging from 563.7 km/s to 155.3 km/s. By fixing the v sin i value to 271.13 km/s and using R = 5.33 R_⊙ (see Table 4), we estimate that the emitting region of the disk expands from 4.9 to 65.0 R_⊙, indicating a significant growth of the disk over time. This growth is likely driven by an enhanced mass loss from the primary star. We computed the Roche lobe radii for the binary system using the formulation of Eggleton (1983) with M₁ = 11 M_⊙ (see Table 4) and M₂ between 0.5 and 0.8 M_⊙ (Tsujimoto et al. 2023; Klement et al. 2024; Huenemoerder et al. 2024) and found that the disk remains well within the Roche radius (R_L ≈ 124–131 R_⊙). This implies that the companion star should not have accreted significant mass from the disk, at least by overflow mass transfer mechanism, as the disk has not reached the Roche lobe boundary (see Fig. 9). Furthermore, considering the orbital separation of 0.96 AU (206.4 R_⊙) proposed by Bjorkman et al. (2002), we conclude that the companion star does not truncate the disk (emitting region of the optical range) of the Be star.

If there is no accretion by overflow mass transfer mechanism, the X-ray emission could come from the interaction between the star and its decretion disk or accretion by other mechanisms. Nazé et al. (2019b) compared the EW Hα and X-ray fluxes from simultaneous data taken from April to December 2018, finding no correlation over time (unusual for a γ Cas star). Later, Huenemoerder et al. (2024) obtained six high-resolution X-ray spectra of π Aquarii with the Chandra/HETG spectrometer from August to October 2022. Their analysis indicates the presence of very hot thermal plasma and coherent variability in the light curve. In five out of six observations, they found a nearly stable hardness ratio (1.50 × 10⁻¹¹ erg/cm² s on August 23 and 27; September 05, 13 and 14), whereas the final observation (October 30) showed a significantly harder (2.54 × 10⁻¹¹ erg/cm² s) and more absorbed spectrum. They attribute the X-ray emission to accretion onto a (magnetic) WD. This means that the magnetic field of the WD could channel material directly from the wind of the Be star or from the edge of its disk without the Roche lobe needing to be full. Moreover, as we concluded before, the disk of π Aquarii is precessing. This means that there are two main scenarios for the accretion to occur: (i) the WD crosses the disk at two points in its orbit, and it is at these moments that it captures material and temporarily increases its X-ray emission, or (ii) if the orbit is very steep, the WD spends most of its time outside the disk. In this case, it can only accrete material from the wind of the Be star, and the accretion rate is much lower than if it had passed through the disk. Following our results, the inclination angle varied ~10° in 20 years, indicating that the first scenario is more probable, and we should expect small increases in X-rays as the WD crosses the disk. On the other hand, as we noted before, FeII 5018 Å covers an emission region larger than Hα, and is the only FeII emission line showing a different shape profile (V/R) compared to the rest of the lines. This could be due to changes in the disk geometry and its alignment with the WD.

Fig. 9 Evolution of the disk radius of π Aquarii derived from Hα emission line measurements using Huang’s law. The disk radius (blue points with error bars) shows significant growth across five epochs, increasing from 4.9 to 64.9 R⊙. The Roche lobe radius (red shaded region, 124–131 R⊙) and orbital separation (green dashed line) indicate that the disk remains well within the Roche lobe and is not truncated by the companion star. Errors represent uncertainties in the Huang radius.

6 Conclusions

A rapid rotator Be star in a binary system, π Aquarii has been studied since the 1930s due to its conspicuous variability. This work aimed to reveal new information on the spectral variability of this star as well as the stellar and disk parameters by characterizing its behavior between several spectral lines that have changed between absorption and emission features over the years. We used observations from such star catalogs as BeSOS and BeSS and obtained the EW and V/R for Balmer, iron, silicon, helium, and other lines. All data are available in an online Excel file.

Our physical parameter estimations based on the spectrum fittings allowed us to obtain good agreement with other works on π Aquarii. Our mass estimates agree with the lower limit of dynamical mass determined by Bjorkman et al. (2002). Moreover, the estimates on temperatures, gravities, and v sin i are compatible with the literature. In this work, we still assumed oblateness of the star due to its fast rotation, obtaining values for polar and equatorial stellar parameters.

Regarding the disk parameters, we studied five epochs that show a large change in the EW in the years 2001, 2011, 2014, 2018, and 2022. We found that the density of the disk decays faster (with n ~ 3.9) in 2001. In this epoch, the disk was not prominent, showing a shell profile in the Hα emission line. Then, in 2011 and 2014, the disk was growing, with a density distribution decaying more slowly (n ~ 3.2). In the last two epochs, 2018 and 2022, the Hα emission line presents a very broad shape with high intensity (almost double the other epochs), and we could not obtain good agreement, mainly because or method does not take into account the electron scattering in the disk and also because we need lower values (n < 3.0) for the n exponent decay of the density distribution (to obtain a much slower decay) to improve our Hα fitting. We assert that the disk of π Aquarii is in a misaligned binary system, and its inclination angles change between ~60° and 70°, going through shell profiles to DPS to triplepeaked to flat-topped profiles, and it now shows a DPS profile (with R≫ V). We expected to implement the SAC method utilized in Arias et al. (2006) to compute the emitting region of the emission lines and compare them with our results using Huang’s relation.

Following our findings, we propose that the (magnetic) WD crosses the disk at two points in its orbit, and it is at these moments it captures material and temporarily increases the X-ray emission. Also, we found that FeII 5018 Å covers an emission region larger than Hα and that it is the only FeII emission line showing a different shape profile (V/R) compared to the rest of the lines, indicating changes in the outer disk probably related to the WD.

The star π Aquarii remains an unsolved case of γ Cas type stars, and further studies are needed to fully understand what is happening with this system. Simultaneous UV and X-ray observations could reveal correlations between changes in the stellar wind and X-ray emissions.

Data availability

The data used in this work are available on Zenodo at https://doi.org/10.5281/zenodo.15039731

Acknowledgements

DC is grateful for the financial support from Fondecyt No. 1231637. CA, MC, and IA are grateful for the financial support from Fondecyt No. 1230131. TS is grateful for the financial support from Fondecyt No. 3230770. MC and CA acknowledge the support from Centro de Astrofísica de Valparaíso. DT acknowledges support from ANID BECAS/DOCTORADO NACIONAL/2024-21240465. This research was partially supported by the super-computing infrastructure of the NLHPC (ECM-02). This project has received funding from the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014–2020) under the Marie Skłodowska-Curie Grant Agreement No. 823734 and is also Co-funded by the European Union (Project 101183150 – OCEANS). This work used BeSOS Catalog, 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 11130702. This work has used the BeSS database, operated at LESIA, Observatoire de Meudon, France: http://basebe.obspm.fr. This work has been possible thanks to the use of AWS-U.Chile-NLHPC credits. Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).

Appendix A Instrument and resolution of BeSS database

Appendix B Corner plots and adjustment of Hα line profiles and SED

Fig. B.1 Corner plot of Hα line π Aquarii observed in 2001. The 2D histograms for each pair of parameters are shown in blue color. In the main diagonal, the PDFs for each parameter are displayed. The colored circles indicate the Spearman coefficient for the correlation between the pairs of parameters. On the upper right inset, the observational data and residuals are plotted. The red line corresponds to the model that maximizes the likelihood, and the blue region represents the reliability interval of 95% constructed from 100 random models sampled by the code.

Fig. B.2 Similar to Fig. B.1 but for Hα observed in 2011.

Fig. B.3 Similar to Fig. B.1 but for Hα observed in 2014.

Fig. B.4 Similar to Fig. B.1 but for Hα observed in 2018.

Fig. B.5 Similar to Fig. B.1 but for Hα observed in 2022.

Fig. B.6 Corner plot of SED π Aquarii. The E(B − V) and parallax were also modeled. The 2D histograms for each pair of parameters are shown in dark orange. In the main diagonal, the PDFs for each parameter are displayed. The colored circles indicate the Spearman coefficient for the correlation between the pairs of parameters. In the upper-right inset, the observational data are in black dots, and residuals are plotted. The dashed red line corresponds to the model that maximizes the likelihood, and the blue region represents the reliability interval of 95% constructed from 100 random models sampled by the code.

Appendix C Evolution of lines in time

Fig. C.1 Comparison of the evolution of Balmer lines Hα, Hβ, and Hγ for the same observation dates.

Fig. C.2 Same as Fig. C.1 but for helium lines.

Fig. C.3 Same as Fig. C.1 but for silicon lines.

Fig. C.4 Same as Fig. C.1 but for iron lines.

References

All Tables

All Figures

	Fig. 1 Dynamical plot of Hα emission line profiles observed between 2001 and 2024.
In the text

Fig. 2 Method used to obtain V/R and DPS for different shape profiles (solid black line). Case 1 (left panel): distinguished DPS profile. The Gaussian fits are represented in pink and red for the violet and red peaks, respectively. The dashed lines indicate the Gaussian fit’s center and the values for both peak intensities. Values for this example are in the legend of the plot. Case 2 (middle panel): the code selects the maximum intensities from each side of the center (dashed vertical lines). Case 3 (right panel): flat emission line. For these cases, the peaks were selected by visual inspection from each side of the geometrical center of the line (dashed vertical lines).

In the text

	Fig. 3 Observed spectrum of the photospheric line HeI 4471 Å (in black). The red line shows the best fit of the ZPEKTR models, and the gray band represents the flux values of the 67 best models.
In the text

	Fig. 4 Equivalent width, V/R, and DPS of Hα plotted against time in MJD. Gray vertical lines denote five specific dates with significant changes in EW, which are discussed in the text. A horizontal line at V/R = 1 is included to highlight the asymmetry of Hα.
In the text

	Fig. 5 Top panel: power spectrum for the V/R presented in units of periods. The date-compensated discrete Fourier transform periodogram is shown as a black solid line, while the Lomb–Scargle periodogram is represented as a blue solid line. Both periodograms have been normalized by their respective means to enhance visibility. Bottom panel: phase of the V/R computed with the period of 84.75 days.
In the text

	Fig. 6 Double-peak separation of iron lines compared to their EW for the penultimate observation date of Fig. C.4. This date corresponds to December 31, 2021. Red vertical dashed lines denote the DPS for the Balmer lines at the same date.
In the text

Fig. 7 Hertzsprung-Russell diagram for π Aquarii. The solid lines represent the evolutionary tracks of rotating stars with solar metallicity and ω = 0.8, as outlined by Georgy et al. (2013) and Granada et al. (2013), the dash magenta line has an initial M = 10.5 M⊙. The filled red circle denotes the findings from this study. The filled blue circle indicates the luminosity and temperature derived from Zharikov et al. (2013), calculated using varying distances. The filled gray circle corresponds to data from Tsujimoto et al. (2023). The filled green circle corresponds to data from Hohle et al. (2010).

In the text

	Fig. 8 Oblateness versus projected velocity of Be stars, as presented by van Belle (2012). The solid circle represents π Aquarii, with data derived from Table 4. Different colors indicate the luminosity class of each star, with the red color representing the star Pi Aquarii.
In the text

Fig. 9 Evolution of the disk radius of π Aquarii derived from Hα emission line measurements using Huang’s law. The disk radius (blue points with error bars) shows significant growth across five epochs, increasing from 4.9 to 64.9 R⊙. The Roche lobe radius (red shaded region, 124–131 R⊙) and orbital separation (green dashed line) indicate that the disk remains well within the Roche lobe and is not truncated by the companion star. Errors represent uncertainties in the Huang radius.

In the text

Fig. B.1 Corner plot of Hα line π Aquarii observed in 2001. The 2D histograms for each pair of parameters are shown in blue color. In the main diagonal, the PDFs for each parameter are displayed. The colored circles indicate the Spearman coefficient for the correlation between the pairs of parameters. On the upper right inset, the observational data and residuals are plotted. The red line corresponds to the model that maximizes the likelihood, and the blue region represents the reliability interval of 95% constructed from 100 random models sampled by the code.

In the text

	Fig. B.2 Similar to Fig. B.1 but for Hα observed in 2011.
In the text

	Fig. B.3 Similar to Fig. B.1 but for Hα observed in 2014.
In the text

	Fig. B.4 Similar to Fig. B.1 but for Hα observed in 2018.
In the text

	Fig. B.5 Similar to Fig. B.1 but for Hα observed in 2022.
In the text

Fig. B.6 Corner plot of SED π Aquarii. The E(B − V) and parallax were also modeled. The 2D histograms for each pair of parameters are shown in dark orange. In the main diagonal, the PDFs for each parameter are displayed. The colored circles indicate the Spearman coefficient for the correlation between the pairs of parameters. In the upper-right inset, the observational data are in black dots, and residuals are plotted. The dashed red line corresponds to the model that maximizes the likelihood, and the blue region represents the reliability interval of 95% constructed from 100 random models sampled by the code.

In the text

	Fig. C.1 Comparison of the evolution of Balmer lines Hα, Hβ, and Hγ for the same observation dates.
In the text

	Fig. C.2 Same as Fig. C.1 but for helium lines.
In the text

	Fig. C.3 Same as Fig. C.1 but for silicon lines.
In the text

	Fig. C.4 Same as Fig. C.1 but for iron lines.
In the text