A study of Be stars in the time domain I. Spectral data and polarimetry

We present the first part of a spectroscopic and polarimetric study on a sample of 58 Be stars that have been measured since 1998. The aim of the study is to understand the timescales of disk variability, formation and dissipation as a function of the properties (mass, luminosity and rotational velocity) of the underlying B star. In this paper we classified the sample based on the presence of emission or absorption of the H 𝛼 line, and the shape of the peak as single or double peak, as well as noting changes between emission and non-emission states. We find a probability of ∼ 0 . 75 percent per year that an object in the sample will undergo such a change. We also present re-derived values of the projected rotational velocities for the sample. When we compare our polarization values with those from the literature, we find that most of the stars do not show a change in the value of the polarization angle, however a small number show significant changes which could be attributed to either disk strength (optical depth) or geometry changes. Finally we show how, by combining the (interstellar corrected) degree of polarization and the projected rotational velocity, we can construct an inclination angle-free parameter that includes the true equatorial velocity. Using this inclination angle-independent parameter we show that the populations of single and double peak stars are indistinguishable, giving further evidence that Be star line profiles are essentially inclination angle driven.


INTRODUCTION
Be stars are, on a fundamental level, B stars that show, or have shown in the past, one or more Balmer lines in emission (Collins 1987).The origin of this emission is currently understood as a circumstellar disk around the star that has been created by material lost by the star (Poeckert & Marlborough 1978;Quirrenbach 1994).
The mechanism that creates this disk is still debated: It is understood that a rapid rotational velocity () plays an important role (Slettebak 1979), although in general they do not rotate at the full break-up velocities (Porter 1996).Mechanisms proposed include single-star radial pulsation (Baade 1988;Rivinius et al. 1998), but it is unclear if this is enough to create the disk (Owocki 2006;Cranmer 2009).Magnetic fields may also play a role (Cassinelli et al. 2002;Brown et al. 2008), but no magnetic fields have been found on Be stars, in contrast with B stars (Wade et al. 2014).Another possibility is that the creation of the disk is due to interaction with a companion star (Kriz & Harmanec 1975;Pols et al. 1991), where the accreting star in the transfer not only increases its mass but also its angular momentum (Hastings et al. 2021).A study from Oudmaĳer & Parr (2010) found that only ∼30% of stars are found in binaries, while a study from Bodensteiner et al. (2020) found that while none of the Be stars studied were in visible binaries with main sequence stars, a large fraction of the fast-rotating B stars were, supporting the idea that binarity may be the cause of the fast rotation and the companion is barely visible anymore after the mass transfer.At the same time, ★ E-mail: a.castanonesteban@2019.ljmu.ac.uk while no main sequence stars have been found as companions of Be systems, subdwarf O and B stars have been found, with characteristics similar to remnant core of a star that did lose most of its envelope (Mourard et al. 2015).A more detailed discussion on the importance of the binarity channel in the creation of the disk around the B star is given by Jones et al. (2022).
In order to help constrain models, observational studies of the properties and variability of the Be phenomenon are likely to be key.For example, understanding the timescales of disk loss and reformation as a property of stellar mass, temperature and equatorial velocity is likely to be important (Barnsley & Steele 2013).Understanding these dependencies is, however, complicated by the effect of inclination.The disk will have a certain inclination, , with respect to the observation plane and the Earth.As a result, the observed rotational velocity is  sin .This value is useful as a lower bound for the rotational velocity, but at the same time the actual rotational velocity  can be much higher.An upper bound on the  will be the critical velocity ( crit ), at which point the rotational and gravitational forces at the equator of the star are balanced.A value often used to see how close to this point a star is the critical fraction  = / crit .
The variability and shape of the emission line has also been a matter of study: It can be a single or double peak (Hanuschik et al. 1988) with variable values of asymmetry between the peaks (Hanuschik et al. 1995) or have the shape of shell lines (Hanuschik 1995).These profiles have been successfully modelled and theoretically reproduced (Silaj et al. 2010).
As well as showing emission line variability, the polarization of Be stars is also variable with time (Coyne & Gehrels 1967).The polarization is caused by radiation scattering due to free electrons (Coyne & Kruszewski 1969) and is a function of the inclination angle.As such the combination of spectroscopic and polarimetric observations can help constrain the true rotational velocity .
In this paper we present the first in a series of new analyses of the properties of a sample of Be stars (Steele et al. 1999) that have been observed spectroscopically for more than two decades.The overall aim of the project is to understand the variability timescales of the objects and how that relates to the properties (mass, luminosity, rotational velocity) of the underlying B star.For this first paper we concentrate on combining new polarimetric observations with the historical spectral data in order to understand the rotational velocity distribution of the sample.This will then be used in future papers to explore the impact on line emission variability at various wavelengths.
In section 2 we describe the properties of the sample and previous work that has been carried out on it.In section 3 we describe the acquisition and data reduction of our spectroscopic and polarimetric dataset and in section 4 make a basic classification of the appearance and variability of the H line profile for each object.Section 5 describes the calculation of projected rotational velocities and critical velocities for the sample.Similarly, section 6 describes the correction for each object of polarization degrees and angles for the effects of interstellar polarization.In Section 7 we combine the results of the previous two sections to derive rotational velocities corrected for the effect of inclination angle.Finally in section 8 we make some concluding remarks.

SAMPLE SELECTION
Our sample is the same as the one introduced in Steele et al. (1999).The 58 stars in the sample were were initially classified as Be stars in the database presented in Jaschek & Egret (1982), and were selected to cover several objects for each spectral and luminosity class for Be stars, to make it as representative as possible.However, because Be stars are variable, some of the stars classified as Be stars in 1982 might not show H emission anymore.As the raw data used in their initial classification is generally unavailable and was compiled from a wide variety of sources, they may also have been misclassified.Several of the stars might also have changed to no emission at some point before or between the observations.Previous work on the sample is described in a series of five papers published between 1999 and 2013.In paper I (Steele et al. 1999), the criteria for the selection of the sample was presented.A basic study of the rotational velocities of the stars was done based on a single spectrum of each object, as well as a reclassification of their spectral type.Paper II (Clark & Steele 2000) studied the same sample plus another 8 sources in the K band (2.05-2.22m).It determined that some objects do not show Br emission ( = 21660 Å), which might indicate no hydrogen emission in general, putting in doubt their classification as a Be star.Different sets of the spectra showed H I, He I, Mg II, Fe II and Na I lines, with He I and Mg II features as a good diagnostic of early spectral type.At the same time, stars that did not show Br emission seemed to have a lower projected rotational velocity.Paper III (Steele & Clark 2001) studied 57 of the stars from paper 1 in the H band (1.53-1.69m).H I Brackett lines were again examined, as well as Fe II.It showed that the analysis of the H band spectra alone only allows for the classification into "early" (B0e-B4e) and "late" (B5e-B9e) types.Paper IV (Howells et al. 2001) studied 52 of the stars from paper 1 with  infrared photometry, separating values for the interstellar reddening and the circumstellar excess, and Table 1.The values given for the defocus of the telescope when using MOP-TOP for the different sources, depending on their apparent magnitude V.

V
Defocus (mm)  > 9.75 0.0 9.75 >  > 6.75 1.0 6.75 >  > 4.9 2.0 4.9 >  > 4.04 2.5 finding a strong correlation between the derived interstellar reddening values and the equivalent width of the interstellar sodium lines, giving confidence in the measured reddening.Paper V (Barnsley & Steele 2013) studied the equivalent width variability of the H line for the measurements taken for 55 stars between 1998 and 2010, finding that stars of earlier types, with higher values of the projected rotational velocity, show a higher degree of variability in the H emission.

Observations and Data Reduction
The SPRAT spectrograph also on the LT (Piascik et al. 2014).The years of observation, as well as the total number of observations, are shown in Table A1.All INT and LT data were re-reduced using the aspired pipeline (Lam et al. 2020) for consistency.
In addition to these spectra, we have polarimetric observations taken with the MOPTOP polarimeter on the Liverpool Telescope (Jermak et al. 2016(Jermak et al. , 2018) ) in October 2022.Our sample is very bright, with apparent magnitudes in the range of  = 4.04 and  = 10.60 (Steele et al. 1999).For this reason, if we try to observe the brightest objects directly, the observation will be instantly saturated.
To avoid this, we defocussed the telescope, with different values depending on how bright the stars are, allowing us to spread the light without saturating the detector.We also used the MOPTOP fast-rotation mode, which completes a cycle of 16 positions in 8 seconds, to avoid saturation as much as possible.The values used for the defocussing are shown in Table 1 to aid future observers of bright objects with the instrument.
Reduction of the MOPTOP data was carried out using aperture photometry routines from astropy (Astropy Collaboration et al. 2022).The extracted photometric counts were converted to Stokes  and  values following the procedure outlined in Shrestha et al. (2020), with corrections for instrumental polarization made using the calibration values presented in Steele et al. (2023).

Additional data
As additional data, we have included the distances and the reddening.The distances have been calculated from parallax values obtained from the Gaia database (Gaia Collaboration et al. 2022), except in 3 cases where they have been taken from the Hipparcos catalog (van Leeuwen 2007) because they were not available or the error of the Gaia values was larger than the actual values.The reddening  (−) was obtained in Howells et al. (2001), where we could find values for 52 of the stars in our data.These values have been added to Table A2.

H𝛼 CLASSIFICATION
The most basic classification that we can make for Be stars from their spectra is taking into account the shape of the H line: some Be stars will show a single peak emission, while others will show a mostly symmetric double peak emission.We show an example of this in Figure 1.At the same time, others will not show any H emission anymore, and will show absorption of the line instead.Even more interesting, with a long time study like this one, we can also see some changing with time.We will discuss this in a quantitative fashion in future papers.However for now we simply classify them based on their overall appearance.In the most complicated cases of the ones with changing emission type, we will choose the most recent observation, as it is the most closely related to their polarization status.
From our sample of 58 stars, 10 never showed any H emission between 1998 and our last observation.16 show a single peak in all of our observations, while 21 show a double peak in all of them.The remaining 11 in our sample have changed with time, with 3 of them showing emission and 8 not showing emission in our last observation.In total, from our last observation, 18 do not show hydrogen emission and 40 show hydrogen emission lines (Figure 2).
Examining the 11 objects that changed with time in more detail, we see that 9 of them at least spent some time as no emission, with another 2 always emitting and changing only between double peak and single peak.We can therefore assume a change rate between emission and absorption of 9/48 in 24 years, or 0.78% of the total sample changing between emission and absorption every year (counting only the objects that we have observed actively changing from our first observation).Similarly if we instead assume all the objects that did not show any absorption in 1998 are changed from their original observations, compiled in Jaschek & Egret (1982) but originally from the 1960s and 1970s, we would have a fraction 19/58 over an estimated time of 45 years (0.73%/year).
Previous determinations of this change rate (McSwain et al. 2008(McSwain et al. , 2009;;Barnsley & Steele 2013;Dimitrov et al. 2018) over shorter time periods (4-12 years) have ranged between 0.3 and 25 %/year.This wide range can be attributed to how the different samples used for the studies were constructed.Using the same sample as here but only between 1998 and 2010 Barnsley & Steele (2013) found 2 changing objects over 55 analyzed Be stars over a period of 12 years, or 0.30%/year.However (McSwain et al. 2008)   The distribution of our sample for the three different types of observed H, in this case, "NP" for "No peak", "SP" for "Single peak", and "DP" for "Double peak", classified from their last observation.
appears that a change rate of no more than a few percent per year is something like the true value.

Calculation of the rotational velocities
One of the most common values for Be stars is the rotational velocity of the central star, .The way that we have chosen to obtain these values is using the He-I absorption lines in our spectra, more precisely the lines at wavelengths 4026 Å, 4143 Å, 4387 Å and 4471 Å.We first fit a Gaussian profile to each of these lines, obtaining this way , the standard deviation of the Gaussian distribution, that then In this case, it is the observation of star HD4180 taken in August 1998.We fit each line to a Gaussian profile and obtained the FWHM from them, and hence the velocity as described in Section 5.The results for each line in this example are vsini( = 4026Å) = 193.7±3.2 km/s, vsini ( = 4143Å) = 190.8±4.9 km/s, vsini ( = 4387Å) = 187.7 ± 4.9 km/s, vsini ( = 4471Å) = 181.9± 3.9 km/s.can be converted into the full width half maximum (FWHM) using: These values of FWHM are not yet a measurement of the rotational velocity itself, since rotation is not the only line broadening mechanism.We therefore convert them using the FWHM- sin  relation derived from observational fits to rotationally broadened line profiles of model atmospheres of a wide sample of early type stars (Slettebak et al. 1975).We will first do a linear fit for the values of FWHM and  of the He-I line ( = 4471Å) shown in the tables 1 and 2 of Slettebak et al. (1975), as shown in Figure 4, and then use this relation to convert our values to rotational velocities, applying a Doppler shift to the FWHM dependent slope to change the values to the other three near He-I lines ( = 4387 Å, = 4143 Å, = 4026 Å).The final equations obtained were: with  () being the FWHM for the line at the corresponding wavelength.
For any particular object the results of the fits to these lines will show certain variation: a number of factors can cause this effect such as the noise in the spectrum or blending with another line.For that reason, we have chosen to calculated the weighted mean and standard deviation of the four lines and reject any values that are more than 4  away from the average value.With this, we have obtained a series of values in time for the rotational velocity of the star.We do not expect any changes in the rotational velocity of a star over the span of just a few decades, so we take the weighted mean of these values for each object.These results are shown in Table 2 and represented  in Figure 5.The red line is a 1:1 line for comparison that allows us to see that our rotational velocities seem to be similar when compared with the lower values of 1999, but lower when compared with the higher values of 1999.The cause of this difference might be that our results are the average of many years of observations, versus only one observation as in 1999.

Calculation of the critical velocity
The critical velocity of a Be star is the theoretical equilibrium velocity at which the central B star would need to rotate over to eject the material on its equator purely by compensating gravity with rotation.We can find the formula for example in Townsend et al. (2004): with M the mass of the central star and   its polar radius.
To obtain these values for our sample, we use classification of the spectral types from Steele et al. (1999) (reproduced in our Table A2) to take stellar luminosities from de Jager & Nieuwenhuĳzen (1987); Cox (2000).Masses and radii were then derived from those luminosity values using the equations shown in Tables 4 and 5 of Eker et al. (2018).
First, we have calculated the masses from the luminosities, following the set of equations: (5) and then we finally calculate the radius of the star using Stefan-Boltzmann law: After obtaining the values of the mass and radius of the star, we can calculate the   and the value of the critical fraction , to show how close the obtained equatorial rotational velocity is to the critical value: We have also calculated the value of  , where in equation 7 we used  instead of .All of these values are listed in Table A2.

Calculation of the polarization degree and angle
After calculating the values for the rotational velocity ( sin ), we move to the calculation of the polarization degree and angle of each band.We first calculate the Stokes parameters following the instructions for the double camera settings as shown in Shrestha et al. (2020).After obtaining the Stokes parameters, we calculated the polarization degree and angle for each band according to the standard formulation.We have derived the polarization degree () and angle () values for the bands ,  and .These results can be found in Tables A3 and A4.As most of the values that we will compare our polarization values (section 7) with are full-range/white-light instead of individual bands, we have calculated the average of the values q and u for the three bands and calculated the average values for the polarization degree ( P) and the angle ( θ) from them.
We note here that MOPTOP polarimetry suffers from a systematic errors of ∼ 0.002 in  and 1 • in  (Shrestha et al. 2020).These systematic errors are not included in the values presented in Tables A3 and A4 but are incorporated into the final interstellar corrected values calculated in the next section.

Separation of the interstellar and intrinsic polarization
We have attempted to separate interstellar polarization from the total polarization.To obtain the values for the interstellar polarization we have first selected stars from the 9286 object catalogue (Heiles 2000) within an angular distance of 5 degrees on the sky from each of our sample stars.We then found the distances to these group of stars using the Gaia catalog.We have separated the values of the Stokes parameters  and  from the catalogue polarization degree and angle for these stars, then performed a linear fit to them to find the values of interstellar  and  with respect the distance.Then we calculated   and   for our sample stars depending on their distance, using the previous linear fit, as shown in Figure 6.The errors for the interstellar polarization are derived from the values of the linear fit.We then subtracted the values of the interstellar  and  for our sample stars from the total values of their Stokes parameters, obtaining the intrinsic values (  and   ), and then calculated the values for both the interstellar and intrinsic polarization degree and angle as shown in Table A7.
In Table 3 we compare our values with the most studied Be stars from our sample, for a historical perspective, while in Table A5 we compare our entire sample with the ones from Yudin (2001).We can see that the polarization degree changes with time, and while the polarization angle is mostly constant (Figure 7), there are a number of objects with significant (> 10 • ) changes over the 70 year time baseline.We highlight these objects in bold in Table 3.We note that large changes in polarization angle are typically associated with large changes in the degree of polarization (either positive or negative).However, we note that the objects that show such angle changes seem to be spread fairly evenly between those which do and do not show changes in the H peak morphology (see last line of table).Overall the changing polarization angle can be interpreted as either a varying ratio of interstellar and circumstellar polarization (which naturally have random polarization angles with respect to each other) as the disk strength varies, or as a change in the disk geometry as a function of time.Such behaviour has, for example, been seen in the interferometric and polarimetric measurements presented by Quirrenbach et al. (1997).We will explore this question on an objectby-object basis in a future paper.

Removal of the inclination dependence
As radiation escapes the Be star and is scattered through the disk, it will get polarized.This polarization has been modeled by Brown & McLean (1977): with P the polarization degree, τ the average electron optical depth,  a shape factor that describes the asymmetry of the disk, and  the inclination.
It is impossible for us to separate these values without knowing the value of the inclination.However, as we have measured previously the rotational velocity  sin , this allows us to obtain instead: This allows us to derive a value for  that is not dependent on the inclination of the star (although there is now an uncertainty related to the mean optical depth instead).Similarly we can multiply the values shown in equation 9 by the critical velocity: obtaining a measurement of the critical fraction .

Testing the population differences
In order to test the validity of the above procedure, we can compare the distribution of the various quantities  sin ,  sin ,  2 / and  2 / Table 3.The values P and θ shown here are the weighted average of their years.For the 1950s, we have used the ones shown in Hall (1958) and Behr (1959), for the 1960s we have used the ones shown in Coyne & Gehrels (1967), Coyne & Kruszewski (1969) and Serkowski (1968), for the 1970 we have used Serkowski (1970), McLean & Clarke (1976), Coyne (1976), Poeckert & Marlborough (1976) and McLean & Brown (1978), for 2001 we have used the ones shown in Yudin (2001), and for 2022 our own calculated values.In bold, the values for the polarization angles that have shown a difference bigger than 10°over the past 70 years.against the H emission classification we made in Section 4. To do this we used Kolmogorov-Smirnov tests to compare the histograms of these quantities for the permutations of no-(NP), single peaked-(SP), and double peaked-emission DP).The results of the tests are shown in Table 4.

Object
In all cases, when considering the distributions that include the sin  inclination term, the distributions for different H emission properties all indicate a different parent population.Similarly the right panel shows the difference between single and double peaks, with double peaks seemingly rotating faster in terms of  sin .However, when we remove the angular dependence on  the inclination ( 2 /) the null hypothesis can not be rejected 9.This can be interpreted as evidence for the idea that the cause of the difference between single and double peak H emission is the inclination of the disk, as predicted by emission line modelling (e.g.Poeckert & Marlborough (1976, 1978), Sigut et al. (2020)).If that theory is taken as accepted, then our method of using polarization to eliminate the dependence on the inclination angle appears to work for the sample distribution.

CONCLUSIONS
Using a sample of 58 stars previously classified as Be, we have used spectra taken over multiple epochs over a 24 year period to classify them based on their their qualitative  line appearance.
We classified these stars in three groups depending on whether they have emission in the shape of a single peak, in the shape of a double peak, or no emission.We found a probability of ∼ 0.75 percent per year that an object would change state between emission and non-emission.We also calculated the projected rotational velocity ( sin ) from multiple measurements of the He-I line profiles.We then compared polarization measurements taken in 2022 with those from previous authors, revealing that the majority of objects show no significant changes in the polarization angle.This is to be expected if this were a mainly geometric parameter based on the system having a stable 3d space orientation.However we do find a small sample of objects which seem to show genuine polarization angle variability which may be interpreted as either varying disk strength (optical depth) or geometry.This will be investigated further in a future paper.

Object
objects in our sample have been observed several times between 1998 and 2022 with multiple instruments.This work includes observations taken in 1998 and 2002 using the IDS spectrograph on the Isaac Newton Telescope (INT), 2009 observations with the medium resolution settings ( ∼ 2500) of the FRODOspec spectrograph on the Liverpool Telescope (LT) (Morales-Rueda et al. 2004), 2010 to 2021 observations with the high resolution ( ∼ 5000) FRODOspec settings, and 2022 observations with the low resolution ( ∼ 350)

Figure 1 .
Figure1.Representative plots of the three types of H peaks in our sample, extracted from the 1998 observations.HD4180 is an example of a single peak, HD22192 is an example of a double peak, and HD164741 is instead an example of absorption in H.The rest wavelength of H is indicated by the vertical line.

Figure 3 .
Figure3.A plot, showing the fit to the helium I lines.In this case, it is the observation of star HD4180 taken in August 1998.We fit each line to a Gaussian profile and obtained the FWHM from them, and hence the velocity as described in Section 5.The results for each line in this example are vsini( = 4026Å) = 193.7±3.2 km/s, vsini ( = 4143Å) = 190.8±4.9 km/s, vsini ( = 4387Å) = 187.7 ± 4.9 km/s, vsini ( = 4471Å) = 181.9± 3.9 km/s.

Figure 6 .
Figure 6.Here we see the polarization coefficient  of nearby stars (in blue) compared to our object, BD+31 4018 (in green).The red line is a fit to the nearby star data.

Figure 8 .
Figure 8. Two histograms showing the distribution of stars for different values and types of emission of the H line with respect to  .

For
example Figure 8 (left) shows the distributions of P and NP objects are different in terms of  sin , with the first showing showing faster  sin  associated with  emission, and the showing slower  sin  associated with  absorption as expected from previous studies.

Figure 9 .
Figure 9.A histogram showing the distribution of stars for single and double peak H line with respect to  2  .

Table 4 .
Results from the KS tests.Distributions which are apparently drawn from the same parent population (i.e.> 0.05) are highlighted in bold.ValuesNP vs P NP vs SP NP vs DP SP vs DP

Table A1 .
The list of observations taken for the respective objects.

Table A6 .
The distance D to the object, and the minimum and maximum distance (d) that we take.The stars included in this range are used for our linear fits to obtain an estimation of the interstellar polarization.The large variation in these ranges is due to two main reasons: loss of linearity at big distances away from the star, or lack of data covering any other suitable range.

Table A7 .
The calculated values for the interstellar and intrinsic polarization degrees and angles.The estimated systematic error from MOPTOP is  ( P) = 0.002 and  ( θ ) = 1 degree.The error of the interstellar polarization values has been obtained from the linear fit to nearby stars, while the error of the intrinsic polarization values has been obtained from the calculations of   and   and the following derivation of   and   .)   (   )    (   )    (   )