The life cycles of Be viscous decretion discs: fundamental disc parameters of 54 SMC Be stars

Abstract

1 INTRODUCTION

In a classical, observational, and quite broad definition, a Be star is a hot, massive star, with a B spectral type (mass ranging roughly from 3 to 17 M_⊙), non-supergiant, whose spectrum has, or had at some time, one or more Balmer lines in emission (Jaschek, Slettebak & Jaschek 1981; Collins 1987). In a more modern and theoretically oriented definition, a Be star is a very rapidly rotating and non-radially pulsating B star that forms a geometrically thin viscous decretion disc (VDD) composed of an outwardly diffusing, viscosity-driven gaseous Keplerian disc that is fed by mass ejected from the central star (Rivinius, Carciofi & Martayan 2013), and a possibly non-negligible line-driven wind (Kee, Owocki & Sundqvist 2016). There is no evidence of large-scale magnetic fields in Be stars (Wade, Grunhut & MiMeS Collaboration 2012; Wade et al. 2016). Fast stellar rotation lowers the effective gravity near the stellar equator and a second mechanism, likely to be stellar pulsation (Rivinius et al. 1998; Baade et al. 2016; Rivinius, Baade & Carciofi 2016), is responsible for pushing this near equatorial matter into orbit in the inner disc. Once in orbit, a viscous mechanism takes place, diffusing matter and angular momentum outwards, thus making the disc grow.

Be stars are usually quite variable in all observables and in several time-scales (days, weeks, months, or even years). The variability observed in Be stars indicates that the injection of matter and angular momentum into the disc is frequently quite erratic, with sudden outbursts of mass injection and periods of no or negligible mass injection. Rivinius et al. (2016) propose a terminology, which will be used here, in which a star that possesses a disc is said to be active and, conversely, when there is no detectable disc, the star is inactive. Two additional terms are used to distinguish the phases of active disc formation (outbursting Be star) and dissipation (dissipating Be star).

It has been demonstrated (e.g. Haubois et al. 2012) that a Be disc fed roughly at a constant rate, and for a sufficiently long time (a few to several years, depending on the value of the viscosity), reaches a quasi-steady state in which the density is nearly constant in time. If the gas temperature is properly taken into consideration, the radial density profile is typically a complicated function of the distance from the star (e.g. Carciofi & Bjorkman 2008). However, a usual approximation is to consider the gas to be isothermal, in which case the density profile assumes a power-law form. This simple steady-state VDD has been successful in describing the main observed features of individual Be discs (Carciofi et al. 2006, 2007; Jones et al. 2008b; Carciofi et al. 2009; Klement et al. 2015, 2017) and samples of Be stars (Silaj et al. 2010; Touhami, Gies & Schaefer 2011; Vieira et al. 2017). However, despite its great success, there are several open and intriguing theoretical questions about the VDD model. Besides the fact that a good description of the mechanism responsible for putting stellar material into orbit is still needed, a complete physical understanding of the forces that drive these discs is also lacking. It has been commonly assumed that the forces operating on the discs are the gravity from the central star and the forces that come from the gradient of pressure and from viscosity, the latter being the one capable of producing torque. Recently, Kee et al. (2016) showed that radiative line forces may also generate a non-negligible torque, at least for gaseous discs with Solar metallicity. In this work, we will proceed with the assumption that line-driven forces are negligible. We mitigate this potential issue by choosing to study Be stars in the SMC, whose low metallicity will greatly decrease the strength of the line forces.

In the alpha-disc formalism, the kinematic viscosity is scaled with the α parameter, defined such that the Rϕ component of the stress tensor is proportional to the gas pressure: W_Rϕ = −αP. The most reliable and direct way of estimating α is to study the time-dependent disc behaviour, where the diffusive effect of viscosity will have clear observational counterparts. Therefore, light curves of long temporal coverage, such as the ones given by microlensing or planetary transit surveys, are excellent instruments to study the dynamical processes in action on the disc, as it builds up and dissipates.

Dynamical studies of Be star viscous discs are still quite scarce. Jones, Sigut & Porter (2008a), using a 1D time-dependent treatment of the alpha-disc and a non-LTE radiative transfer code, studied the temperature and density profiles of a dynamical disc and their respective H α line profiles. Haubois et al. (2012) studied the theoretical photometric effects of time variable mass injection rates on the structure of the disc also using a 1D time-dependent treatment of the alpha-disc, associated with the Monte Carlo radiative transfer code hdust (Carciofi & Bjorkman 2006, 2008). Carciofi et al. (2012), by fitting these dynamical models to a dissipating portion in the light curve of the Be star 28 CMa (which passed from an outbursting phase, that lasted from 2001 to 2003, to a dissipating phase at the end of 2003), estimated the value of the α parameter for the Be disc of 28 CMa to be α = 1.0 ± 0.2. Later, however, it was realized that a proper consideration of the previous history of the disc must be taken into consideration even when fitting the dissipating portion of the light curve, and this quite high value has been revisited to be closer to α = 0.2 (Ghoreyshi & Carciofi 2017).

Another intriguing issue regarding Be stars is how they acquired such high rotation rates (typically 80 per cent of breakup, Rivinius et al. 2013). As the rotating B star evolves, core contraction and internal angular momentum redistribution generally tends to enhance surface angular rotation (Ekström et al. 2008; Granada et al. 2013). Another scenario (e.g. Pols et al. 1991) would involve a past mass-transfer phase in a binary system, during which the primary donates mass and angular momentum to the secondary. The left-over of such a system would be a fast-spinning Be star (the former mass gainer) and a subdwarf O or B star (sdO/sdB, the former mass donor). Regardless of how they were spun-up, it has been proposed (Krtička, Owocki & Meynet 2011) that the discs of Be stars may provide natural mechanisms for removing large quantities of angular momentum from the fast rotating stars, preventing them to reach the rotation critical limit. The evolutionary models of Granada et al. (2013) assumed the appearance of completely formed viscous discs every time their models reached a near-critical rotation. The mass density and the rate of angular momentum loss of their discs were roughly similar to the ones estimated by Vieira et al. (2017), who modelled the spectral energy distribution (SED) of 80 Be stars using the VDD model, provided that values of α of at least a few tenths were assumed in both approaches.

The main objective of this paper is to build upon the previous dynamical studies of Be discs and provide, for the first time, a detailed study of the temporal evolution of a large sample of Be stars from the Small Magellanic Cloud (SMC). By studying a large sample of stars, we may begin to answer several open questions related to Be stars and their discs, namely: (i) What is the typical value of the viscosity in these discs? (ii) Is there any significant evidence for a dependence of α with parameters such as the density of the disc, the spectral type of the star, etc.? (iii) What are the typical rates of mass and angular momentum loss in these stars? To reach these goals we developed a new method for modelling the light curves of Be stars, described in Sections 2 and 3. The sample of studied light curves is described in Section 4, and the model results are discussed in Section 5, followed by the conclusions.

2 VISCOUS DECRETION DISCS AROUND BE STARS

The optical light curves of early-type Be stars (with spectral type ranging roughly from B0 to B4) are usually quite variable in time-scales of days to years, with amplitudes of up to tenths of a magnitude (Rivinius et al. 2013). The majority of them show very irregular variability. Most present clear single bump-like features, characterized by a fast rise in brightness followed by a slower fading. Frequently, between the brightening and the fading phases some sort of plateau of nearly constant brightness is seen. Sometimes these bumps are reversed, so that an initial fast fading is followed by a slow recovering of the stellar brightness. Two examples of light curves showing these dips are presented in Fig. 1, taken from OGLE-II (Udalski, Kubiak & Szymanski 1997) and OGLE-III (Udalski et al. 2008) data, for two Be star candidates from the SMC, based on the selection made by Mennickent et al. (2002). Object SMC_SC1 75701 shows two bumps, while SMC_SC6 128831 shows a dip.

Figure 1.

Two light curves, in photometric bands V (green) and I (red), selected from the OGLE-II and OGLE-III photometric surveys. Above: Light curve of SMC_SC1 75701, showing two bumps. Below: Light curve of SMC_SC6 128831, showing a dip. The pair of vertical dotted straight lines near JD-2450000 = 2000 separates OGLE-II from OGLE-III data. The measurements shown in purple are assumed to represent the inactive (discless) brightness level of the Be star. Their mean is given by the horizontal purple straight lines. The pairs of vertical orange straight lines bracket our visually selected bumps. These bumps are modelled in Figs 8 and 9.

Open in new tabDownload slide

These bumps and dips resemble the photometric features shown in Haubois et al. (2012, e.g. their fig. 14), where a circumstellar viscous disc builds up as a result of a mass injection into the disc at a constant rate, and then dissipates after the mass injection ceases. Haubois et al. (2012) studied several disc feeding scenarios (constant, cyclic, and outburst) and their photometric counterparts, and demonstrated that the bumps are disc formation/dissipation events of active Be stars seen at near pole-on inclination angles (⁠|$i\lesssim 70\deg$|⁠), while the dips are associated with near edge-on Be stars (often called shell stars, |$i\gtrsim 70\deg$|⁠). The inclination angle plays an important role on how the stellar brightness is modified by the presence of disc because, in the second case (edge-on), the disc is seen projected against the stellar disc, thus causing an attenuation of the stellar radiation. This attenuation does not happen for the pole-on case, where the net effect of the disc is to increase the optical brightness as a result of free–bound and free–free radiation from the gas (Gehrz, Hackwell & Jones 1974; Vieira, Carciofi & Bjorkman 2015).

The orbital velocity of the disc, |$v$|_ϕ, is assumed to be Keplerian (⁠|$v_\phi = v_K = v_\mathrm{orb}\tilde{R}^{-{1}/{2}}$|⁠, with |$v$|_orb = (GM/R_eq)^1/2 and |$\tilde{R}=R/R_\mathrm{eq}$|⁠, where R_eq is the stellar equatorial radius). This assumption holds as long as the gravitational force is much larger than the force due to the pressure gradient, which is generally true for distances from the star of dozens to a few hundreds of stellar radii (Okazaki 2001), as long as (c_s/|$v$|_K)² ≪ 1.

In addition, mass can flow away from the disc through its boundaries. Mass can fall back into the star through the inner boundary at R_eq, or it can leave the system at an outer boundary, R_out. The outer boundary can be interpreted as the limiting radius of the disc due to a binary companion (e.g. Okazaki et al. 2002) or due to the photoevaporation of the disc (e.g. Okazaki 2001). We consider that all mass that eventually reaches the stellar equator R_eq is totally absorbed. The same is assumed for the mass that eventually reaches an outer radius R_out (see below for the definition of R_out). Therefore, the boundaries consist of Σ(R_eq, t) = Σ(R_out, t) = 0.

The solution of equation (5) is scaled in time by the time-scale parameter, which controls how fast matter is distributed throughout the disc and, consequently, its observational counterparts. It follows that by fitting observed light curves of Be stars, this parameter can be estimated, and once other parameters are known (R_eq, M, T_disc), the α value can be determined (e.g. Carciofi et al. 2012).

Before moving on to the modelling of bump-like events such as the ones of Fig. 1, we introduce in the next section some important parameters of the dynamical discs implied in the equations shown above.

2.1 Dynamical disc parameters

We may extend equation (11) to the general case of a time-dependent |$\dot{M}_\mathrm{inj}(\tilde{\tau })$|⁠, which would define, by the same equation, a time-dependent asymptotic surface density, |$\Sigma _0(\tilde{\tau })$|⁠. The function |$\Sigma _0(\tilde{\tau })$| therefore is just another way of specifying the history of mass injection from the star into the disc. It has, however, the advantage of being a surface density, which is a quantity that may be determined from, e.g. SED analyses, in contrast to the mass injection rate and the radius of mass injection, which are parameters that cannot be observationally determined.

The steady-state solution (equation 10) shows that, in the wider domain |$\tilde{R}_\mathrm{inj} \le \tilde{R}\le \tilde{R}_\mathrm{out}$|⁠, the density profile of the disc is not altered if |$\tilde{R}_\mathrm{inj}$| is changed, provided that |$\dot{M}_\mathrm{inj}$| is also changed in order to maintain Σ₀ fixed, according to equation (11). In fact, we verified that the time-dependent solutions of equation (5) in the domain |$\tilde{R}_\mathrm{inj} \le \tilde{R}\le \tilde{R}_\mathrm{out}$| are negligibly affected by the particular choice of |$\tilde{R}_\mathrm{inj}$| or |$\dot{M}_\mathrm{inj}(\tilde{\tau })$|⁠, as long as the quantity |$\Sigma _0(\tilde{\tau })$| is kept fixed. This is a consequence of the fact that the dynamical solutions reach a near steady state very quickly in the vicinity of the injection radius (Haubois et al. 2012). Furthermore, provided that mass is injected not too far from the stellar photosphere (i.e. assuming |$\tilde{R}_\mathrm{inj} \gtrapprox 1$|⁠), the domain |$1\le \tilde{R}\le \tilde{R}_\mathrm{inj}$| is much narrower than the region where the continuum visual flux of Be stars is generated (Carciofi 2011), which means that the emission from this region can be ignored. Consequently, we conclude that |$\Sigma _0(\tilde{\tau })$| (with the assumption that |$\tilde{R}_\mathrm{inj}\gtrapprox 1$|⁠) is a much better parameter for describing the mass injection history of the disc than the pair of parameters |$\dot{M}_\mathrm{inj}(\tilde{\tau })$| and |$\tilde{R}_\mathrm{inj}$|⁠.

The time-dependent solutions of equation (5) generally show that, for Be stars dynamically feeding the disc but still far from steady state, the mass flux close to |$\tilde{R}_\mathrm{inj}$| has absolute values of the order of |$\left(-\mathrm{\partial} M/\mathrm{\partial} t\right)_\mathrm{typ}$|⁠, defined by equation (11). Therefore, we refer to this quantity as the typical decretion rate, which depends on parameters relatively easy to estimate from SEDs of Be stars.

In our simulations, since the values of |$\dot{M}_\mathrm{inj}$| and |$\tilde{R}_\mathrm{inj}$| are of no interest, and the value of |$\tilde{R}_\mathrm{out}$| is quite uncertain, we arbitrarily chose |$\tilde{R}_\mathrm{inj}=1.017$| and |$\tilde{R}_\mathrm{out}=1000$| (we discuss below how this choice of |$\tilde{R}_\mathrm{out}$| might affect our results). Equation (11) therefore shows that the typical decretion rate is much smaller than the mass injection rate |$\dot{M}_\mathrm{inj}(\tilde{\tau })$|⁠. In our case, the typical decretion rate is only |$8.46\times 10^{-3} \dot{M}_\mathrm{inj}$|⁠. This means that the majority of the injected mass flows inwards and is absorbed by the inner boundary at the stellar equator, and only a small remaining fraction of the injected mass is responsible for the growth of the disc. These results were first obtained from SPH simulations of Be discs by Okazaki et al. (2002), who found that only about 0.1 per cent of the injected material flows outward, as a direct result of their choice for R_inj.

2.2 The mass reservoir effect

It is important to stress that the solution |$\Sigma (\tilde{R}, \tilde{\tau })$| is shaped not just by the mass injection rate |$\Sigma _0(\tilde{\tau })$| at the specific instant |$\tilde{\tau }$|⁠, but by the whole mass injection history before the instant |$\tilde{\tau }$|⁠. Therefore, the advantage of studying the relatively isolated bumps like the ones exemplified in Fig. 1, which started after a clear inactive phase, is that there is no disc present when the bump starts developing; thus, no previous history of mass injection has to be taken into account in the beginning of the modelling.

The light curves of several Be stars show that the duration of the build-up phase, which we refer to as the build-up time, is variable between Be stars and even between different bumps from the same star, ranging from a few days to years. The following phase of disc dissipation, however, contrary to the build-up phase, depends of the previous history of mass injection. For this reason, the modelling of the dissipation phase must not be disconnected from the modelling of the build-up phase that happened before it.

One of the main consequences of this fact is the mass reservoir effect (see also Ghoreyshi & Carciofi 2017). Basically, discs that had a longer build-up phase necessarily transported more matter and angular momentum outwards and created a larger external reservoir of mass and angular momentum in its outer regions, which usually extend far beyond the first few stellar radii where the visible photometric observables are formed. It is common, for instance, that some bumps reach plateaus during the build-up phase. The plateau indicates that the density in the inner disc has reached near-steady-state values and, consequently, there is little photometric variation in the visible wavelengths. The outer disc, however, will likely be far from steady state and thus will continue to increase in density and mass. When mass injection ceases, the dissipation phase begins. Re-accretion occurs and, due to the more massive outer disc, the inner disc remains relatively denser for a longer time. This makes the dissipation of the disc appear slower in the observed light curves. Conversely, a disc that had a small build-up time would dissipate much faster.

The importance of the mass reservoir effect can be assessed by the reevaluation of the α parameter in 28 CMa by Ghoreyshi & Carciofi (2017). Carciofi et al. (2012) modelled the 2003 dissipation phase of 28 CMa by considering a very long previous build-up time, and found that a high value of α was necessary (1.0 ± 0.2) to match the observed dissipation rate. Ghoreyshi & Carciofi have shown that when the previous build-up phase is properly accounted for in the modelling, the value of α required to match the dissipation rate is much smaller (0.21 ± 0.05).

3 A MODEL GRID OF DISC FORMATION AND DISSIPATION EVENTS

In this section, we describe the method we developed for fitting the light curves associated with events of disc formation and dissipation. The method consists of pre-computing a large grid of dynamical models of the time-dependent disc structure, covering the entire range of observed scenarios (Section 3.1) and performing the radiative transfer in these models to produce synthetic light curves (Section 3.2). The observed light curves are then fitted by the synthetic one using the procedure described in Sections 3.3 and 3.4.

3.1 Dynamical model grid

For building a comprehensive grid of dynamical models that are solutions of equation (5), we used the definitions of Section 2 that allow us to write the solution |$\Sigma (\tilde{R}, \tilde{\tau })$| in terms of |$\Sigma _0(\tilde{\tau })$| and the dimensionless parameter |$\tilde{\tau }$|⁠.

As discussed in Section 2.2, the advantage of studying relatively isolated bumps like the ones exemplified in Fig. 1, which started after a clear inactive phase, is that there is no previous history of mass injection to be taken into account for the modelling, so that during build-up the shape of the curve is controlled solely by |$\Sigma _0(\tilde{\tau })$| and |$\tilde{\tau }$|⁠, while for dissipation the previous disc build-up time should also be considered (Section 2.2). By using the time parameter (equation 7) instead of the physical time, our dynamical models are independent of the specific physical parameters M, T_eff, R_eq, and α(t) of the Be star under consideration (equation 6). Also, from the linearity of equation (5), it follows that multiplying |$\Sigma _0(\tilde{\tau })$| by some constant results in the solution |$\Sigma (\tilde{R}, \tilde{\tau })$| multiplied by the same constant. Consequently, only one value of Σ₀ during the build-up phase is necessary.

For our grid of dynamical models, we therefore assume that our Be stars start discless. At instant |$\tilde{\tau }=0$|⁠, mass injection into the disc begins at an arbitrary constant rate (Σ₀ > 0) that lasts until |$\tilde{\tau }=\tilde{\tau }_\mathrm{bu}$|⁠, which we refer to as the scaled build-up time, since it is related to the above-mentioned build-up time, but scaled by the time-scale parameter. After that (⁠|$\tilde{\tau }>\tilde{\tau }_\mathrm{bu}$|⁠), mass injection no longer occurs (Σ₀ = 0) and the disc dissipates. In Appendix A, we further discuss the properties of these dynamical models.

We chose 11 values of |$\tilde{\tau }_\mathrm{bu}$|⁠, listed in Table 1. Since the time-scale parameter (equation 6) is roughly given by ∼(100–200)/α d for early Be stars in the main sequence with α ≲ 1, these values correspond to real build-up times of at least 15 d, which brackets the observed build-up times of the sample described below (Section 4, Table 5). In this study, we decided not to model the bumps with observed build-up times lower than about 15 d, usually referred to as flickers (Keller et al. 2002).

3.2 Radiative transfer models

Having selected a set of suitable hydrodynamic bump models, the next step is to produce photometric light curves of these models. The radiative transfer part of the problem requires a stellar model, which will be the primary source of radiation. The stellar model depends on the physical parameters M, R_eq, and T_eff, which were left unspecified in the dynamical model grid. In addition, three other parameters must be specified: the viewing angle, i (i = 0 means pole-on orientation), the distance to the star, d, and the interstellar reddening.

One important feature of the central stars of Be stars is that they are fast rotators. Fast rotation causes the star to be oblate, with hotter poles and colder equatorial regions. Rotation is specified by the ratio of the rotation velocity at the equator to the Keplerian velocity at the equator, W = |$v$|_eq/|$v$|_orb. The ratio between the equatorial radius to the polar radius is given by R_eq/R_pole = 1 + W²/2 for a Roche-shaped star. All these parameters evolve in time as a consequence of stellar evolution, and Be stars can be found in luminosity classes from V to III (Rivinius et al. 2013). We adopt the Geneva evolutionary tracks (Georgy et al. 2013) to determine R_eq and T_eff given M and the age in the main sequence, t_MS.

The current version of hdust allows for a spheroidal rotationally oblate star, with the latitude-dependent surface temperature being given by |$T_\mathrm{surf}\propto g_\mathrm{eff}^\beta$| (Carciofi et al. 2008). Here, the coefficient β(W) is calculated by fitting a straight line to the gradient |$\mathrm{\partial} \ln T_\mathrm{surf}/\mathrm{\partial} \ln g_\mathrm{eff}$| given by the flux theory of Espinosa Lara & Rieutord (2011). For the disc scale height (equation 9), we assume an isothermal disc with T_disc = 0.6T_eff, where T_eff is the effective temperature of the star, defined by T_eff = (L_*)^1/4(σS_*)^−1/4, with S_* being the surface area of the star.

In order to generate synthetic absolute magnitudes from the computed SEDs, we used the standard BVRI Johnson–Cousins passbands from Bessell (1990) and the Vega flux from Castelli & Kurucz (1994) as standard of calibration.

A grid of model light curves was computed using the 11 dynamical models described in Section 3.1, with 8 different values of Σ₀ (third column of Table 1). For each of these disc models, radiative transfer models were calculated with hdust at 17 different time parameters (not shown in the table) and 15 equally spaced values of cos i (second column). This whole process was done for three different stellar models (‘Star 1’, ‘Star 2’, and ‘Star 3’, first column of Table 1, according to the stellar models of Georgy et al. 2013). Details on the stellar models are given in Table 2. They were chosen to represent early B-type stars from the SMC (Z = 0.002), in the middle of their life in the main sequence, with the rotation parameter given by the mean value obtained for Be stars (W = 0.81, Rivinius, Štefl & Baade 2006). In the sixth column of Table 2, we present the values of ατ (equation 6) for the discs of these stars (with the assumption that T_disc = 0.6T_eff). In short, a single light curve is specified by taking one element of each column of Table 1. The end result was a grid of 3 × 15 × 8 × 11 = 3960 light curves, for each of the BVRI bands.

A grid of inactive (discless) stellar models was also calculated. Because these models can be computed much faster than the bump models, we were able to cover a much finer grid of stellar parameters (Table 3), aiming at a better determination of the stellar parameters. The grid is composed by models for 13 different masses (second column), 5 different rotation rates (third column), 6 equally spaced values for the age in the main sequence (forth column), and 10 equally spaced values of cos i (fifth column), resulting in a total of 13 × 5 × 6 × 10 = 3900 photometric models for each of the BVRI bands.

3.3 Empirical law

In order to facilitate the comparison of the synthetic light curves (Section 3.2) with the observed ones (Section 4), we developed two empirical laws that match quite closely the synthetic light curves for build-up and dissipation. The usefulness of these formulae will become clear in the next section.

In our discussion of the features of the light curves, it is useful to separate them in three groups: (i) pole-on light curves, of stars seen at small inclination angles (⁠|$0\le i \lesssim 70\deg$|⁠), which should statistically correspond to the majority of the observed light curves; (ii) edge-on light curves, of shell stars (⁠|$i \approx 90\deg$|⁠); and (iii) intermediate light curves, of stars seen at intermediate angles (⁠|$70 \lesssim i \lesssim 85\deg$| – the extension of this intermediate region varies depending the photometric band under consideration and will be defined below). Pole-on light curves show an increase in apparent brightness, due to the additional flux coming from the disc. Conversely, edge-on light curves show a decrease in apparent brightness, due to obscuration of the star by the disc. The intermediate case shows the smallest variations in apparent brightness, and frequently the light curve has a more complicated shape, as it is influenced by variable amounts of disc emission/absorption.

|$\Delta X_\mathrm{d}(\tilde{\tau })/\Delta X_\mathrm{d}^0$| is a function that goes from 1 to 0 as |$\tilde{\tau }-\tilde{\tau }_\mathrm{bu}$| goes from 0 to ∞. In Appendix B, we show examples of light curves that accompany the conclusions drawn on this section.

The values of |$\Delta I_\mathrm{bu}^\infty$| for our grid are shown in Fig. 2, plotted against cos i. The values for the BVR bands show qualitatively similar patterns to the ones presented in this figure. Each panel shows the results for a different star, and all 11 values of Σ₀ (Table 1) are represented in the figure by different colours and symbols. The curves monotonically increase with cos i, starting with negative values at edge-on orientation and reaching a maximum for pole-on viewing. The angle for which |$\Delta X_\mathrm{bu}^\infty = 0$|⁠, where the disc excess emission is exactly matched by the absorption of photospheric light by the disc, depends both on the density scale (as shown in the figure) and (most importantly) on the band pass.

Figure 2.

The values of |$\Delta I_\mathrm{bu}^\infty$| versus cos i for our grid. From top to bottom, the results are for Star 1, Star 2, and Star 3, respectively. Purple, blue, green, orange, and red circles correspond to Σ₀ = 2.50, 1.85, 1.37, 1.01, 0.75 g cm⁻². Purple, blue, and green triangles correspond to Σ₀ = 0.56, 0.41, 0.30 g cm⁻². Vertical dotted lines define the region of intermediate angles for the I band (⁠|$73 \lesssim i \lesssim 84\deg$|⁠).

Open in new tabDownload slide

An analysis of our model grid allowed us to determine the ranges in inclination angle for which the light curve displays the intermediate behaviour described above. They were determined by visual inspection of our model grid, as the angles for which the light curves present more complex shapes (see e.g. the I-band light curves seen at |$i=81.8\deg$| and |$i=77.6\deg$| in Fig. B1 of Appendix B). Their adopted values are 53–|$78\deg$|⁠, 60–|$78\deg$|⁠, 66–|$84\deg$|⁠, and 73–|$84\deg$| for the BVRI bands, respectively.

Fig. 2 also shows that the excesses increase a little when moving from a low- to a high-mass star, for discs with the same other features. This is a consequence of the fact that the stellar flux relative to the disc flux increases with the luminosity of the star.

Figure 3.

Comparison of the empirical law (black straight lines) with the computed light curves for build-up and dissipation. Top panels: log –log derivative of |$1-\Delta I_\mathrm{bu}(\tilde{\tau })/\Delta I_\mathrm{bu}^\infty$| versus |$\Delta I_\mathrm{bu}(\tilde{\tau })/\Delta I_\mathrm{bu}^\infty$|⁠. Bottom panels: log –log derivative of |$\Delta {I}_\mathrm{d}/\Delta {I}_\mathrm{d}^0$| versus |$1-\Delta {I}_\mathrm{d}/\Delta {I}_\mathrm{d}^0$|⁠, for values of |$\tilde{\tau }_\mathrm{bu}$| equal to 0.45, 1.5, 6, and 30. The results are shown for Star 2 at two inclination angles: |$i=0\,\deg$| (left), and |$i=90\,\deg$| (right). The green, red, and blue curves correspond to Σ₀ equal to 1.37, 0.75, and 0.41 g cm⁻².

Open in new tabDownload slide

Representative I-band values of ξ_bu and ξ_d for our grid are shown in Figs 4–6 (the values of ξ_bu and ξ_d for the BVR bands show qualitatively similar patterns to the ones presented in these figures). Each row shows the results for a different star. The left (right) panels are for edge-on (pole-on) models. The values of ξ are directly related to the rate of photometric variations: the smaller the ξ, the slower the variation (see equations 17 and 18).

Figure 4.

The I-band values of ξ_bu for our grid versus cos i. Left: Edge-on models. Right: Pole-on models. From top to bottom, the results for Star 1, 2, and 3, respectively. The markers are the same as in Fig. 2.

Open in new tabDownload slide

Figure 5.

Selected I-band values of ξ_d for our grid versus cos i. The scaled build-up time was fixed to |$\tilde{\tau }_\mathrm{bu}=2.25$|⁠. Left: Edge-on models. Right: Pole-on models. From top to bottom, the results for Star 1, 2, and 3, respectively. The markers are the same as in Fig. 2.

Open in new tabDownload slide

Figure 6.

Selected I-band values of ξ_d for our grid versus |$\tilde{\tau }_\mathrm{bu}$|⁠. Left: Edge-on models with |$i=90\,\deg$|⁠. Right: Pole-on models with |$i=0\,\deg$|⁠. From top to bottom, the results for Star 1, 2, and 3, respectively. The markers are the same as in Fig. 2.

Open in new tabDownload slide

Fig. 4 exemplifies the strong variation of ξ_bu with i, for the pole-on case. This is probably due to the fact that, in the build-up process, the density grows from the inside out (see Appendix A), which means that the optically thick part of the disc (an expanding pseudo-photosphere) forms first, and the optically thin part takes longer to be built. Since the optical excess of the disc is given by ΔX ∝ F_thickcos i + F_thin, it follows that, as we move from pole-on to edge-on angles, only the optically thick contribution (proportional to cos i) varies. As a consequence, the observed rate of increase in flux moves from being more to less optically thick dominated.

In the dissipation process, the density rapidly adjusts to a self-similar dissipation pattern in the inner disc (Appendix A). Therefore, the flux from the disc is the result of the decrease and disappearance of the optically thick region – transformed into an optically thin region – and the diminishing of the whole optically thin region. The pole-on values of ξ_d (right panels of Fig. 5) are affected by cos i to a less extent, when compared to the values of ξ_bu. In the dissipation process, by the same reasoning applied to the build-up process, since the optically thick emission is attenuated by the effect of cos i, its disappearance has a reduced effect for more inclined discs and therefore the disc should apparently disappear at a slower rate.

The values of ξ_d also show great variation with the asymptotic surface density (Figs 5 and 6). More specifically, increasing Σ₀ results in a light curve with a slower decay rate. This is probably due to different levels of saturation in the optically thick region. The denser the optically thick region, the bigger its optical depth and the greater the amount of time for it to turn into an optically thin region.

In addition, Fig. 6 shows that ξ_d strongly depends on the scaled build-up time. As expected from the mass reservoir effect (Section 2.2), increasing |$\tilde{\tau }_\mathrm{bu}$| results in smaller values of ξ_d, which implies slower decay rates.

From the above an important conclusion can be drawn: viscosity is not the only parameter affecting the rate of photometric variations in a Be light curve. The stellar parameters, the asymptotic surface density, as well as the inclination angle, all affect the observed shape of the light curve. Thus, extracting α from light curves, one of the main goals of this paper, cannot be done without some knowledge about these parameters.

3.4 Fitting pipeline

So far, our model light curves were given in terms of the adimensional time parameter, |$\tilde{\tau }$|⁠. Thus, an equation is necessary to transform from the physical time t to |$\tilde{\tau }$|⁠, in order to connect the real light curves to our simulated ones.

Our goal is to fit an observed light curve with equation (21), in order to obtain, in a self-consistent way, all the stellar and disc parameters of interest. For that, the following chain of procedures is adopted:

• Find a light curve of a Be star that contains at least one clear inactive phase and one complete photometric bump.

• Obtain the magnitudes X_* at the inactive phase. Subtract these magnitudes from the light curve and obtain the excesses ΔX(t).

  Without a clear inactive phase, it is not possible to obtain the pure photospheric level (e.g. the horizontal purple straight lines in Fig. 1) and, consequently, it is not possible to know how much of the observed bumps represent the disc contribution to the total flux. In addition, the photometric bump must contain a completely identified build-up phase, from which the instants t₁ and t₂ can be extracted, and a considerable extension of the dissipation phase.

• Fit equation (21) to the selected bumps, obtaining the coefficients |$\Delta X_\mathrm{bu}^\infty$|⁠, C_bu and C_d, as well as the times t₁ and t₂ for the onsets of build-up and dissipation.

• Transform the magnitudes at the inactive phase, X_*, to absolute magnitudes, M_X*, by correcting for the distance to the star and reddening at each observed band.

  Given the theoretical dependence of the coefficients in equation (21) on the stellar parameters, the absolute magnitudes are required to estimate the stellar parameters (M, W, and t/t_MS). From them, the parameter ατ (see equation 6) can be estimated. Clearly, if the stellar parameters are known from some other way (e.g. by spectroscopic analysis), this requirement is no longer necessary. Unfortunately, this is not the case for our sample.

• Estimate the stellar parameters, the geometric parameter (cos i) and the bump parameters (Σ₀, α_bu, and α_d, for each bump) that best reproduce the fitted stellar (M_X*) and bump (⁠|$\Delta X_\mathrm{bu}^\infty$|⁠, C_bu, and C_d) parameters [see equations (22) and (23) and the parameters |$\Delta X_\mathrm{bu}^\infty$|⁠, ξ_bu, and ξ_d from Section 3.3].

In practice, the above process involves several complications (e.g. estimating the goodness of the fit) that are described in the next section.

3.5 Fitting using a MCMC sampling

The task of fitting the measured stellar absolute magnitude (M_X*) and bump parameters (⁠|$\Delta X_\mathrm{bu}^\infty$|⁠, C_bu, and C_d) for estimating the model parameters – step 5, above – was done using the Markov Chain Monte Carlo (MCMC) sampling technique. We used the python MCMC sampler emcee (Foreman-Mackey et al. 2013). The code samples a large collection of models by varying all model parameters within a pre-specified range. The sampler provides a distribution of model parameters according to a posterior distribution p(model∣data) ∝ L(data∣model)π(model), where L and π are the likelihood and the prior distributions, respectively.

In our fitting procedure, there are 4 + 5N_bumps model parameters for each light curve containing N_bumps identified bumps. There are three stellar parameters (M, t/t_MS, W) and one geometric parameter (cos i), and, for each bump in the light curve, there are five parameters: the initial times of the build-up and dissipation phases (t₁ and t₂), the asymptotic surface density (Σ₀), and the viscosity parameters during the build-up and dissipation phases (α_bu and α_d).

For the parameter sampling, we have chosen hundreds of ‘walkers’,¹ proportional to the number of 4 + 5N_bumps model parameters. For a randomly chosen set of parameters sampled by emcee in the course of the simulation, the corresponding stellar and bump observables are calculated by a multidimensional linear interpolation of the model grid. During the raffle of parameters the prior probability was set to zero if one of the values were sampled outside of the allowed range of a given parameter (Tables 1 and 3). The simulation consists of two steps, the so-called burn-in phase and the sampling phase. We verified that 1000 iterations in the burn-in phase were sufficient for the convergence of all our models. For each parameter, the best-fitting values were chosen to be the median of distribution of the posterior probabilities, with upper (lower) uncertainties estimated from the differences between |$84\,\,\rm{per\,\,cent}$|(⁠|$16\,\,\rm{per\,\,cent}$|⁠) of the sample and the median

In order to test our fitting routine, we applied it to synthetic light curves with levels of astrophysical noise and uncertainties similar to the observed light curves used in this work. In general, the parameters used to generate the synthetic light curves were fairly recovered, with the exception of cos i and α_bu, for which a strong correlation is expected on theoretical grounds (Fig. 4). As expected, when multiband data was used the errors in the derived parameters were smaller. This test indicates that Σ₀ and α_d can be reliably estimated from the light curves, while cos i and α_bu less so.

In Section 4, we present a selection of light curves of Be stars from the SMC, and measure their stellar and bump quantities (steps 1 to 4 of the pipeline). Later, in Section 5, we apply step 5 of the pipeline, as described in this section, in order to estimate the relevant parameters of the selected Be stars.

4 OGLE LIGHT CURVES OF BE STAR CANDIDATES

Mennickent et al. (2002) selected roughly one thousand Be star candidates from the SMC, by studying light-curve variations using the OGLE-II data base (Udalski et al. 1997). They classified the morphologies found in the light curves into four categories. The majority of the light curves (⁠|${\sim } 65\,\,\rm{per\,\,cent}$|⁠) belonged to their type-4 category, composed by the light curves showing irregular and non-periodic variations. These light curves should correspond to Be stars showing episodes of mass injection more complicated than the simple build-up followed by dissipation scenario described in Section 3.

Most interesting for us is the type-1 group (⁠|${\sim }13\,\,\rm{per\,\,cent}$| of the sample), composed by light curves that show single sharp or hump-like bumps, like the bumps of the light curve of SMC_SC1 75701 (Fig. 1). These bumps should be the result of single nearly continuous episodes of mass injection followed by dissipation of the disc, like the theoretical scenario explored in Section 3.

The type-2 group (⁠|${\sim }14\,\,\rm{per\,\,cent}$| of the sample), containing light curves showing high and low plateaus, also has some interesting cases for our purposes. High plateaus are usually the photometric result of a longer build-up process in which a near steady state has been reached in the inner disc. The low plateaus are frequently the portions of the light curve during inactive phases.

Sabogal et al. (2005) selected roughly two thousand Be star candidates from the LMC and classified their light curves into the same four categories described by Mennickent et al. (2002). Previously, Keller et al. (2002) also studied light curves from the LMC using the MACHO survey. They spectroscopically analysed a subsample of their Be star candidates and found that |$90\,\,\rm{per\,\,cent}$| of them were Be stars. They also classified morphologically their light curves in a slightly different manner. Their so called bumper events and ‘flicker events’ more or less correspond to the bumps of type-1 light curves, but also to features of the more irregular type-4 light curves. The bumpers have duration of a few hundred days, while the flicker events are faster, with durations of a few dozens of days. Dips like the one exemplified by the light curve of SMC_SC6 128831 (Fig. 1) were called ‘fading events’. The frequency of these events was quite smaller than the bumpers, in accordance to the picture that fading events are associated with the less numerous shell stars.

Paul et al. (2012) studied the spectral properties of stars from the catalogues of Mennickent et al. (2002) and Sabogal et al. (2005). For the candidates from the SMC, they found that the majority of type-1 and type-2 light curves belong to early B-type stars with emission features characteristic of circumstellar material (Paul et al. 2012). Therefore, these light curves are very likely to be from Be stars.

In this work, we selected light curves from the catalogue of Be star candidates from the SMC of Mennickent et al. (2002). In order to have light curves of a longer time baseline, we combined OGLE-II data with OGLE-III (Udalski et al. 2008). Due to a calibration issue between OGLE-II and OGLE-III, namely a shift in the zero-points present in some of the light curves, it was necessary to find inactivity intervals in both the OGLE-II and OGLE-III portions of these light curves to measure and correct the problem.

The light curves were visually inspected according to the criteria of item 1, Section 3.4, i.e. light curves with at least one clear inactive phase and one bump. In this initial work we focused on well-behaved light curves with clear bumps. We also avoided the short events (flickers, with build-up times ≲15 d), due to the fact that most of them are poorly sampled. The end result was a sample of 54 stars, containing 81 selected bumps, shown in Table 5. In the table, horizontal lines separate the data for each of the 54 stars. Each row in the table contain the data for each of the 81 selected bumps. The fifth and sixth columns in the table contain the beginning and ending of the selected inactive interval for the light curve. The seventh, eighth, and ninth columns contain the B_*V_*I_* magnitudes obtained at the inactive phase for the light curve. Due to the nature of the OGLE survey, the B_*V_* are not available for all sources. The 11th column contains the bands that were considered in the fitting process of the specific bumps, depending on the availability of measurements in each band. The last two columns are initial visual estimates of t₁ and t₂, which were used as input for emcee.