# Critical study of the distribution of rotational velocities of Be stars

J. Zorec, Y. Frémat, A. Domiciano de Souza, F. Royer, L. Cidale, A.-M. Hubert, T. Semaan, C. Martayan, Y. Cochetti, M. Arias, et al.

► To cite this version:

J. Zorec, Y. Frémat, A. Domiciano de Souza, F. Royer, L. Cidale, et al. Critical study of the distribution of rotational velocities of Be stars. Astronomy and Astrophysics - A&A, 2016, 595, pp.A132. 10.1051/0004-6361/201628760. hal-02425734

HAL Id: hal-02425734

https://hal.science/hal-02425734

Submitted on 14 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Critical study of the distribution of rotational velocities of Be stars

## I. Deconvolution methods, effects due to gravity darkening, macroturbulence, and binarity*

J. Zorec<fnref n="1" />,<fnref n="2" />, Y. Frémat<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, A. Domiciano de Souza<fnref n="4" />, F. Royer<fnref n="5" />, L. Cidale<fnref n="6" />,<fnref n="7" />, A.-M. Hubert<fnref n="5" />, T. Semaan<fnref n="8" />, C. Martayan<fnref n="9" />, Y. R. Cochetti<fnref n="6" />,<fnref n="7" />, M. L. Arias<fnref n="6" />,<fnref n="7" />, Y. Aidelman<fnref n="6" />,<fnref n="7" />, and P. Stee<fnref n="4" />

Received 21 April 2016 / Accepted 8 June 2016

## ABSTRACT

*Context.* Among intermediate-mass and massive stars, Be stars are the fastest rotators in the main sequence (MS) and, as such, these stars are a cornerstone to validate models of structure and evolution of rotating stars. Several phenomena, however, induce under- or overestimations either of their apparent $V\sin i$, or true velocity $V$.

*Aims.* In the present contribution we aim at obtaining distributions of true rotational velocities corrected for systematic effects induced by the rapid rotation itself, macroturbulent velocities, and binarity.

*Methods.* We study a set of 233 Be stars by assuming they have inclination angles distributed at random. We critically discuss the methods of Cranmer and Lucy-Richardson, which enable us to transform a distribution of projected velocities into another distribution of true rotational velocities, where the gravitational darkening effect on the $V\sin i$ parameter is considered in different ways. We conclude that iterative algorithm by Lucy-Richardson responds at best to the purposes of the present work, but it requires a thorough determination of the stellar fundamental parameters.

*Results.* We conclude that once the mode of ratios of the true velocities of Be stars attains the value $V/V_c \approx 0.77$ in the main-sequence (MS) evolutionary phase, it remains unchanged up to the end of the MS lifespan. The statistical corrections found on the distribution of ratios $V/V_c$ for overestimations of $V\sin i$, due to macroturbulent motions and binarity, produce a shift of this distribution toward lower values of $V/V_c$ when Be stars in all MS evolutionary stages are considered together. The mode of the final distribution obtained is at $V/V_c \approx 0.65$. This distribution has a nearly symmetric distribution and shows that the Be phenomenon is characterized by a wide range of true velocity ratios $0.3 \le V/V_c \le 0.95$. It thus suggests that the probability that Be stars are critical rotators is extremely low.

*Conclusions.* The corrections attempted in the present work represent an initial step to infer indications about the nature of the Be-star surface rotation that will be studied in the second paper of this series.

**Key words.** stars: emission-line, Be – stars: rotation

## 1. Introduction

A classical Be star is a non-supergiant B-type star whose spectrum has, or had shown at some time one or more Balmer lines in emission (Jaschek et al. 1981; Collins 1987). This definition selects objects in the main-sequence evolutionary phase (MS) underscoring thus the “active” character of the Be phenomenon. These stars have the faculty of creating their own circumstellar environment in which emission lines are produced. This is

in contrast to early-type supergiants in which emission in the Balmer lines is a “passive” consequence of their extended atmosphere. Even when some classic Be stars are believed to survive as such in the bright giant phase (Negueruela 2004), this definition excludes the heterogeneous group of objects presenting the B[e] phenomenon, which likely manifests at different stages of stellar evolution (Lamers et al. 1998). Be stars have long been known as rapid rotators (Struve 1931) and, accordingly, these stars are assumed to create their circumstellar disk or envelope thanks mainly to this rapid rotation. General properties of Be stars and mechanisms to build up their circumstellar envelopes (hereafter CE) or disks are reviewed by Rivinius et al. (2013).

<footnote n="1">¹ Sorbonne Universités, UPMC Université Paris 6 et CNRS UMR 7095, Institut d'Astrophysique de Paris, 75014 Paris, France
e-mail: zorec@iap.fr</footnote>

<footnote n="2">² CNRS UMR 7095, Institut d'Astrophysique de Paris, 98bis bd Arago, 75014 Paris, France</footnote>

<footnote n="3">³ Royal Observatory of Belgium, 3 avenue Circulaire, 1180 Bruxelles, Belgium</footnote>

<footnote n="4">⁴ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS UMR 7293, Lagrange, 28 avenue Valrose, 06108 Nice Cedex 2, France</footnote>

<footnote n="5">⁵ GEPI, Observatoire de Paris, PSL Research University, CNRS UMR 8111, Université Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92190 Meudon, France</footnote>

<footnote n="6">⁶ Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque S/N, 1900 La Plata, Argentina</footnote>

<footnote n="7">⁷ Instituto de Astrofísica La Plata, CONICET, 1900 La Plata, Argentina</footnote>

<footnote n="8">⁸ Geneva Observatory, University of Geneva, chemin des Maillettes 51, 1290 Sauverny, Switzerland</footnote>

<footnote n="9">⁹ European Organization for Astronomical Research in the Southern Hemisphere, Alonso de Cordova 3107, Vitacura, Santiago de Chile, Chile</footnote>

<footnote>* Full Tables 1 and 4 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/595/A132</footnote>
The rapid decrease of the source function of Balmer lines in two domains of effective temperatures, those larger than $T_{\text{eff}} \approx 27\,000$ K, and lower than $T_{\text{eff}} \approx 12\,000$ K, makes the corresponding emission signatures produced by CE or disks too weak to be detected in the respective spectral type domains (Zorec et al. 2007). Adding to this the steep drop of the IMF for masses $M \gtrsim 15\,M_{\odot}$, the observed frequency of Be stars displays two maxima and becomes significant only in a limited range of stellar masses, i.e., $3 \le M/M_{\odot} \le 30$. Nothing precludes, however, that continuous and variable winds implying mass-loss rates $\dot{M} \sim 10^{-11}-10^{-9}M_{\odot}/\text{yr}$ (Snow 1981) and episodic mass ejec-tions up to $\dot{M} \sim 10^{-7}M_{\odot}/\text{yr}$ that are typical in Be stars (Brown & Wood 1992; Hanuscihk et al. 1993; Floquet et al. 2000; Hubert et al. 2000; Keller et al. 2002; Mennickent et al. 2002; de Wit et al. 2006; Meilland et al. 2006) be present in more- and less-massive rapidly rotating stars than in those of B spectral type considered in the present paper. The conclusions drawn in this work may also concern stars in these two mass extremes.

Regarding the critical rotation, statistical analyses (Cranmer 2005; Frémat et al. 2005) conclude that Be stars are under-critical rotators. Studies of individual stars, based on nonradial pulsation modes (Cameron et al. 2008; Saio 2013) and rotational rates inferred from interferometric data (Meilland et al. 2012; Domiciano de Souza et al. 2014), produce estimates of $\Omega/\Omega_c$ ranging from 0.60 to 0.98, which correspond to linear velocities ratios $0.44 \le V/V_c \le 0.88$ ($V_c$ is the critical equatorial linear velocity) or to ratios of centrifugal to gravity accelerations $0.14 \le a_c/a_g \le 0.70$. These rates suggest that the surface rotational velocities of Be stars are undoubtedly very high, but they are not necessarily critical. It is then important to inquire how rapid the rotation of Be stars actually is and how much complementary energy is then necessary to release the mass required to build up CEs or disks. On the other hand, these stars represent a paradigm for rapid rotators and, as such, they are the cornerstone to validate the physical input in models of the structure and evolution of rotating stars.

In most cases, the stellar rotation is apprehended through the $V\sin i$ parameter. Nevertheless, a number of phenomena act either to underestimate or increase the $V\sin i$ of Be stars systematically. In the present paper we only pay attention to the following phenomena: (1) the gravity darkening effect (hereafter GD), which may carry a systematic underestimation of the $V\sin i$ (Stoeckley 1968; Townsend et al. 2004; Cranmer 2005; Frémat et al. 2005); (2) more or less chaotic macroscopic velocity fields present in atmospheres of active stars, called macroturbulence, and since Unsöld & Struve (1949) and Huang & Struve (1953) are known to contribute a nonrotational broadening component of spectral lines; (3) under certain conditions, the orbital motion increases the apparent value of $V\sin i$ because this phenomenon can concern a large fraction of Be stars that are members of binary systems (Kriz & Harmanec 1975; Harmanec 1987; Gies 2000; Chini et al. 2012; Nasseri et al. 2013; Oudmaijer & Parr 2010).

The final distribution $\Phi(u)$ of velocity ratios $u = V/V_c$ obtained in the present work will be analyzed in the following paper (Paper II, Zorec et al. 2016) to determine the degree of surface differential rotation that may affect the $V\sin i$ parameters and is consistent with the function $\Phi(u)$ derived here.

The present work (Paper I) is organized as follows. In Sect. 2 presents the stellar sample and the data used in the paper. The Sect. 3 is devoted to the discussion of pros and cons of the most commonly used methods to derive the distribution of rotational velocities corrected for the inclination angle factor $\sin i$. The effect of the macroturbulence and the presence of binary

components on the distribution of velocity ratios are studied in Sects. 4 and 5, respectively. A discussion and comments on the results obtained are summarized in Sect. 6.

## 2. Data

### 2.1. The studied stellar sample

The stellar sample used in this work has 233 Galactic classical Be stars, which have been previously studied by Chauville et al. (2001), Frémat et al. (2005), Zorec et al. (2005), Frémat et al. (2006) and Levenhagen & Leister (2006). They are listed in Table 1. This sample was chosen not only because the $V\sin i$ parameters were determined using very similar methods, but also because the spectra were corrected for CE emission perturbations.

As the discussion in Sect. 3 requires that all stars had fundamental parameters determined, a significant part of this section is dedicated to their estimation. Details concerning the acquisition and reduction of data proper are found in the papers cited above.

### 2.2. Fundamental parameters

Rapidly rotating stars undergo geometrical deformations and the concomitant GD. They can then be characterized with two types of fundamental parameters:

1. Apparent fundamental parameters. Following Frémat et al. (2005), the directly observed parameters are called apparent fundamental parameters, which describe only the average properties of the photosphere in the hemisphere projected toward the observer. In the present work we deal with four independent observed apparent quantities:

(i) $T_{\text{eff}}^{\text{app}}$, is the apparent effective temperature, which corresponds to the effective temperature of a classical nonrotating, plane-parallel model of stellar atmospheres without GD that describes the observed spectrum in a limited wavelength domain.

(ii) As for the effective temperature, $\log g_{\text{eff}}^{\text{app}}$ is the apparent effective surface gravity of the observed stellar hemisphere determined with classical plane-parallel model atmospheres.

(iii) $L_{\text{app}}$, is the apparent bolometric luminosity emitted by the aspect-angle dependent stellar hemisphere.

(iv) $V\sin i$, is the apparent rotational velocity of a rigidly rotating star that accounts for the observed spectral line broadening in the visual spectral range when GD is not taken into account.

2. Parent nonrotating counterpart parameters. To make the calculation of quantities, such as the stellar mass, stellar age, or the inclination angle $i$, easier, in Sect. 3.5.1 we determine the stellar parent nonrotating counterpart parameters. These are called, in short, the *pnrc fundamental parameters* that correspond to the associated nonrotating stars having the same mass as the studied objects (Frémat et al. 2005).

Although the parameters ($T_{\text{eff}}$, $\log g$, $\log L/L_{\odot}$, $V\sin i$) of our stellar sample were discussed several times (see Chauville et al. 2001; Frémat et al. 2005, 2006; Zorec et al. 2005), since then new data have been obtained for some of them. These data enabled us to redetermine their bolometric luminosities. We briefly detail the characteristics of the apparent parameters adopted in the present work.
### 2.2.1. The apparent effective temperature and surface gravity

The apparent effective temperature $T_{\text{eff}}^{\text{app}}$ and surface gravity $\log g_{\text{eff}}^{\text{app}}$ of stars employed in this work were determined in two different ways: (a) from spectrophotometric observations in the 3500 $\le \lambda \le 4600$ Å wavelength interval; and (b) using high resolution spectra in the blue spectral range. In Table 1 gives the sources of the observed parameters, wavelength intervals observed, characteristics of the spectra, and instruments used.

a) The stars in Chauville et al. (2001), Frémat et al. (2005) and Zorec et al. (2005) have apparent effective temperatures and gravities determined from the $(\lambda_1, D)$ parameters defined in the BCD (Barbier, Chalonge, Divan) spectrophotometric system (Chalonge & Divan 1952). For a number of these stars, we added new BCD observations obtained by Cochetti (2014), which are indicated with “4” in Table 1. We redetermined the $T_{\text{eff}}^{\text{app}}$ for these Be stars using the calibrations given in Zorec et al. (2009), while for $\log g_{\text{eff}}^{\text{app}}$ we used those listed in Zorec (1986) and Divan & Zorec (1982a).

The $(\lambda_1, D)$ parameters describe the energy distribution around the Balmer discontinuity (hereafter BD), which is roughly from $\lambda 3500$ Å to $\lambda 4500$ Å. Explanations on the characteristics, use and advantages of the BCD $(\lambda_1, D)$ parameters to determine $T_{\text{eff}}^{\text{app}}$ and $\log g_{\text{eff}}^{\text{app}}$ of Be stars can be found in Zorec & Briot (1991), Cidale et al. (2001) and Gkouvelis et al. (2016). To indicate average MK spectral types, we employ the $S_{70}$ parameter (Zorec et al. 1983) that in the BCD $(\lambda_1, D)$-classification diagram (Chalonge & Divan 1952) corresponds to the value of the BD determined by the intersection of the $\lambda_1 - 3700 = 70$ Å line with the curve of intrinsic color gradient $\Phi_{rb}$ passing through the stellar $(\lambda_1, D)$ point in the BCD $(\lambda_1, D)$-classification diagram (Chalonge & Divan 1952). The $S_{70}$ identifies curves of constant effective temperatures (Zorec et al. 2009) and nicely correlates with the spectral classification parameter “s” defined by de Jager & Nieuwenhuijzen (1987).

b) The $T_{\text{eff}}^{\text{app}}$ and $\log g_{\text{eff}}^{\text{app}}$ from Frémat et al. (2006) and Levenhagen & Leister (2006) were determined by fitting spectra in the blue spectral range and/or ratios of H, He, Si, and N lines in the 3900 $\le \lambda \le 5000$ Å spectral range. Depending on the intensity of the emission in the Balmer lines, corrections to the spectra were carried out for veiling effect. Details on this technique can be found in Ballereau et al. (1995) and Semaan et al. (2013). For those stars where effective temperatures and gravities were determined using both spectrophotometric and spectroscopic data, we adopted the weighted averages, where the weights were assigned using the individual uncertainties of parameters determined in both methods.

A statement detailing the uncertainties of the fundamental parameters introduced by the use of the BCD quantities $(\lambda_1, D)$ can be found in Zorec & Briot (1991) and Zorec et al. (2009). A comparison of effective temperatures determined with the BCD $(\lambda_1, D)$ parameters for B stars without emission with those obtained in the literature using other methods is shown in Appendix A.1.

### 2.2.2. The apparent bolometric luminosity

We call $L_{\lambda}(\eta, i)$ the monochromatic luminosity emitted at the wavelength $\lambda$ by the hemisphere seen under the inclination

angle $i$ of a star rotating at a velocity characterized by the ratio $\eta$ defined as

$$ \eta = \frac{a_c}{a_g} = \left( \frac{\Omega_c}{\Omega_g} \right)^2 \left[ \frac{R_c(\eta)}{R_c} \right]^3, \quad (1) $$

where $a_c$ is the centrifugal acceleration and $a_g$ the gravitational acceleration at the equator; $\Omega_c$ and $\Omega_g$ are the actual and critical stellar angular velocities at the equator, respectively; and $R_c(\eta)$ and $R_c$ are the rotationally modified stellar equatorial and critical radii, respectively.

According to the data we have at our disposal, we write the apparent bolometric luminosity $L^{\text{app}}$ as

$$ L^{\text{app}} = \int_{0}^{\lambda_{\text{UV}}} L_{\lambda} \, d\lambda + \int_{\lambda_{\text{UV}}}^{\lambda_{\text{IR}}} L_{\lambda} \, d\lambda + \int_{\lambda_{\text{IR}}}^{\infty} L_{\lambda} \, d\lambda, \quad (2) $$

where $\lambda_{\text{UV}}$ is the shortest wavelength at which was observed the energy distribution of a star in the far-UV; $\lambda_{\text{IR}}$ is the longest wavelength observed in the far-infrared (IR). The parameters $\lambda_{\text{UV}}$ and $\lambda_{\text{IR}}$ are different from star to star, according to the satellites or instruments that recorded the data. To make easier the estimate of $L^{\text{app}}$, Eq. (2) is rewritten as

$$ L^{\text{app}} = [1 + \Delta_{\text{L}}(\eta, i, M, t)] \int_{\lambda_{\text{UV}}}^{\lambda_{\text{IR}}} L_{\lambda} \, d\lambda \quad (3) $$

$$ \Delta_{\text{L}}(\eta, i, M, t) = \Sigma_{i} \int_{\Delta_{\text{L}}\lambda_{i}} F_{\lambda} \, d\lambda / \int_{\Delta_{\text{L}}\lambda_{i}}^{\lambda_{\text{IR}}} F_{\lambda} \, d\lambda, $$

where $M$ is the stellar mass and $t$ its age. The spectral regions that are not attainable by direct observation are completed with $\Delta_{\text{L}}(\eta, i, M, t)$; these are determined with model-atmosphere fluxes $F_{\lambda}$ of rotating stars calculated with FASTROT (Frémat et al. 2005). The completion factor $\Delta_{\text{L}}$ in Eq. (3) is adjusted consistently, according to the iterated parameters $(\eta, i, M, t)$, as explained in Sect. 3.5.1.

In $L^{\text{app}}$, only the factor $\int_{\lambda_{\text{UV}}}^{\lambda_{\text{IR}}} L_{\lambda} \, d\lambda$ is a genuine observed quantity, which was determined following the same procedure as described in detail by Zorec et al. (2009). Although the use of this method is only justified for rotating stars without emission, it can also be applied to objects with emission lines, such as Be stars, provided that some modifications are introduced that take the characteristics of their energy distributions into account.

The energy distribution in the visual spectral range of Be stars is variable, displaying flux excesses or deficiencies according to the sub-wavelength domain and/or emission phase B normal, Be proper, or Be shell (Divan et al. 1982, 1983; Divan & Zorec 1982b; Zorec & Briot 1991; Moujtahid et al. 1998, 1999). To extract the nonvariable photospheric component $L_{\lambda}^*$ from the luminosity $L_{\lambda}$ emitted by the star+CE system, we have written

$$ L_{\lambda} = L_{\lambda}^{*}(1 + \epsilon_{\lambda}), \quad (4) $$

where $\epsilon_{\lambda}$ represents the contribution from the CE, so that the apparent luminosity $L_{\text{app}}$ becomes

$$ L^{\text{app}} = [1 + \Delta_{\text{L}}(\eta, i, M, t)] \int_{\lambda_{\text{UV}}}^{\lambda_{\text{IR}}} \frac{L_{\lambda}}{1 + \epsilon_{\lambda}} \, d\lambda. \quad (5) $$

For $\epsilon_{\lambda}$ we used the expression obtained in Moujtahid (1998) (Moujtahid et al. 1999, 2000a,b), which accounts for the energy distribution of Be stars from $\lambda \sim 3000$ Å to the near-IR and is a function of the following parameters: $r = R_{*}/R_{\text{env}}$ is the inverse of the radius of the effective layer representing the circumstellar environment; $B_{\lambda}(T_{\text{env}})$ is the source function of this layer given by the Planck function at a temperature $T_{\text{env}}$; $\tau_{\lambda} = \tau_{e} + \tau_{v}(\lambda/\lambda_{v})^3$
<table><caption>**Table 1.** Studied sample of Be stars and their apparent fundamental parameters.</caption><thead><tr><th rowspan="2">HD</th><th rowspan="2">T<sub>eff</sub><sup>app</sup> ± σ<sub>T<sub>eff</sub></sub></th><th>log</th><th>log</th><th rowspan="2">E(B-V) ± σ<sub>E</sub></th><th rowspan="2">(V sin i)<sub>app</sub> ± σ<sub>Vsin i</sub></th><th rowspan="2">Notes</th></tr><tr><th>g<sub>eff</sub><sup>app</sup> ± σ<sub>log g</sub></th><th>L<sup>app</sup>/L<sub>☉</sub> ± σ<sub>log L</sub></th></tr><tr><th></th><th></th><th>K</th><th>dex</th><th></th><th>mag</th><th></th><th></th></tr></thead><tbody><tr><td>144</td><td>12 180 ± 170</td><td>3.80 ± 0.10</td><td>2.964 ± 0.026</td><td>0.08 ± 0.03</td><td>125 ± 6</td><td>1, f</td></tr><tr><td>4180</td><td>15 520 ± 360</td><td>3.54 ± 0.31</td><td>3.426 ± 0.037</td><td>0.11 ± 0.02</td><td>195 ± 10</td><td>1, f</td></tr><tr><td>5394</td><td>35 210 ± 2130</td><td>4.18 ± 0.31</td><td>5.190 ± 0.055</td><td>0.13 ± 0.07</td><td>432 ± 34</td><td>1, f</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>330 950</td><td>29 870 ± 3390</td><td>4.45 ± 0.38</td><td>2.578 ± 0.120</td><td>0.69 ± 0.13</td><td>60 ± 4</td><td>2, 4, d</td></tr></tbody></table>

**Notes.** 1: stars studied in Chauville et al. (2001), Frémat et al. (2005), Zorec et al. (2005); 2: stars studied in Levenhagen & Leister (2006); 3: stars studied in Frémat et al. (2006); 4: stars with new BCD observations (Cochetti 2014); f: stars with integrated fluxes; and d: stars with spectroscopic distances. The complete version of this table is available at the CDS.

(Moujtahid et al. 1999) is the optical depth at $\lambda$ of the equivalent layer, where $\tau_e$ is the electron scattering component, and $\tau_V$ accounts for the bf+ff opacity at $\lambda_V \approx 5500$ Å. These quantities are only used as fitting parameters so that they are not a specific physical interpretation in any case. For those stars for which we had enough data, the CE parameters were determined in the same way as described in Moujtahid et al. (2000a,b) using the Moujtahid et al. (1998) catalog that lists changes with time of: (a) the total BD $D = D_* + d$ of Be stars, where $D_*$ is the photospheric BD and $d \gtrsim 0$ is the second BD due to the CE; and (b) the flux and color excess at $\lambda_V \approx 5500$ Å. The explicit form of $\epsilon_\lambda$ used and the parameters ($r$, $T_{env}$, $\tau_e$, $\tau_V$) obtained for the ten program Be stars with rather strong emission characteristics are given in Appendix B.

The monochromatic fluxes used to determine $L_\lambda$ and calculate the integral in Eq. (5) are from the spectrophotometric data obtained in the far-UV by the TD1, ANS satellites (Jamar et al. 1976; Macau-Hercot et al. 1978; Wesselius et al. 1982) and from the IUE satellite low resolution spectra (CDS compilation). In the visible and near-IR domains the data are from the BCD spectrophotometric data collection, the 13-color photometry (Johnson & Mitchell 1975) and various UVB, uvby, UVBRIJHK, or JHK photometric catalogs (CDS database). The far-IR data are from the ISO infrared sources (CDS). In Table 1 we identified with “f” those stars where $log L^{app}/L_{\odot}$ was estimated by integrating the fluxes from $\lambda_{UV}$ to $\lambda_{IR}$.

For some program stars, we did not find the required far-UV and/or far-IR spectrophotometry. In these cases we could not do better than conform with the apparent bolometric luminosity estimated from the bolometric absolute magnitude $M_{bol} = M_V + BC(T_{eff})$ with $M_{bol}^{\circ} = +4.742$ mag, where $M_V$ is the visual absolute magnitude and $BC(T_{eff})$ the bolometric correction (Flower 1996; Torres 2010; Nieva 2013). The magnitude $M_V$ was determined using the apparent magnitude V of an emissionless phase.

No matter which method is used, to estimate $log L^{app}/L_{\odot}$ we need to know the ISM extinction $A_\lambda = k_\lambda E(B-V)$, where $k_\lambda$ is the ISM extinction law given by Cardelli et al. (1988, 1989) and $E(B-V)$ is the ISM determined in several ways: a) from a smoothed distributions of the absorption $E(B-V)$ against the distance of stars without emission located in as narrow as possible intervals of equatorial spherical coordinates ($\Delta\alpha, \Delta\delta \le 1^\circ$) around each studied Be star [see examples of its use in Aidelman et al. (2012)]; b) by rectifying the energy distributions near the $\lambda$ 2200 Å absorption bump (Beeckmans & Hubert-Delplace 1980; Briot & Zorec 1987) using TD1, ANS far-UV spectrophotometric observations, and/or IUE low resolution spectra, or through the fitting parameters of the 2200 UV bump given by

Guerlter et al. (1982), Friedemann et al. (1983), Friedemann & Roeder (1987) as used in Zorec & Briot (1991). Table 1 lists the adopted average values of all independent determinations of $E(B-V)$ for each program Be star and their respective $1\sigma$ dispersion.

The transformation from the apparent bolometric fluxes to bolometric luminosities is finally realized using HIPPARCOS parallaxes (van Leeuwen 2007). For stars where the trigonometric parallaxes are negative or nonexistent, we had to estimate spectroscopic distances. To this end, we used the diagrams of the long-term variability of the apparent magnitude V made by Moujtahid et al. (1998), which enabled us to seize the value corresponding to Be star emissionless-like phase. The stars with spectroscopic distances are identified with “d” in Table 1.

Appendix A.2 gives a short statement on the systematic deviations that may characterize our estimates of bolometric luminosities. We use for this average MK spectral-luminosity classes of emissionless B-type stars.

### 2.2.3. The apparent rotational velocity

A projected rotational velocity $V\sin i$ determined as detailed below is considered here as an apparent quantity because it can be more or less underestimated due to the GD (Stoeckley 1968; Townsend et al. 2004; Frémat et al. 2005).

For those program stars where $T_{eff}^{app}$ and $log g_{eff}^{app}$ were redefined, as mentioned in Sect. 2.2.1, we estimated the $V\sin i$ parameter using the robust CERN-MINUIT algorithm to fit the spectra (minimization package available at CERN<fnref n="1" />). Otherwise, the $V\sin i$ parameters were determined by Chauville et al. (2001), Frémat et al. (2006), and Levenhagen & Leister (2006), who took similar precautions and approaches to determine these parameters. The main characteristics of the adopted $V\sin i$ parameters

a) are based on the interpretation of He I 4471 and Mg II 4481 lines corrected for veiling effect;

b) give a measurement of the line broadening as produced by rigidly rotating atmospheres without GD;

c) are determined by comparing the observed line profiles with synthetic spectra to avoid uncertainties produced by the use of a unique limb darkening coefficient over the wavelengths inside the spectral lines (Collins & Truax 1995).

When several determinations existed for a given star, we adopted the weighted average. These averages and their $1\sigma$ deviations are given in Table 1.

<footnote n="1">¹ http://wwwasdoc.web.cern.ch/wwwasdoc/minuit/minmain.html</footnote>
## 2.3. The critical linear velocity

Gathering Be stars, and rapid rotators in general, into spectral type groups and luminosity classes makes little sense because both characteristics are aspect-angle dependent and do not clearly reflect either the true stellar mass or its evolutionary stage. In order to keep a statistically reliable set of stars and minimize the bf uncertainties introduced by differences in masses and/or evolutionary stages on the studied velocity distributions, the calculations performed in this work use the $V\sin i$ parameters normalized to the critical linear rotational velocity $V_c$. We adopted the following definition for the critical linear rotational velocity $V_c$:

$$V_c = 436.7 \left[ \frac{M/M_{\odot}}{R_c(M, t)/R_{\odot}} \right]^{1/2} \text{ km s}^{-1}, \quad (6)$$

whose use for the studied Be stars is justified in Appendix C. The stellar masses needed to estimate $V_c$ were interpolated in the evolutionary tracks calculated by Ekström et al. (2012) with metallicity $Z = 0.014$ suited for the solar environment (Asplund et al. 2009). The entries to these models are the *pnrc* effective temperatures and bolometric luminosities defined in Sects. 2.2.1 and 2.2.2. They also determine the *pnrc* stellar radius $R_0(M, t)/R_{\odot}$ used to obtain $R_c(M, t)/R_{\odot}$, where the ratio $R_c(M, t)/R_0(M, t)$ comes from two-dimensional (2D) models of rotating stars (Zorec et al. 2011; Zorec & Royer 2012).

The mass and radius of the parent nonrotating object needed to estimate $V_c$ cannot be estimated in a simple way from the apparent effective temperature and apparent bolometric luminosity because in rapidly rotating stars both are deeply marred by rotational effects. Moreover, they are both more or less aspect-angle dependent. On the other hand, the apparent $V\sin i$ parameters have to be corrected for underestimation induced by the rapid rotation of stars. Since this correction, as well as its use to pass from a distribution of apparent $V\sin i$ to the corresponding one of true velocities $V$, both depend on the *pnrc* parameters, we first define in Sect. 3 these corrections and the required deconvolution methods. We postpone the formal determination of the *pnrc* quantities to Sect. 3.5.1 where, in particular, we consistently obtain the mass and radius of the parent nonrotating object needed to estimate $V_c$.

We can wonder whether systematic differences exist between our estimates of $V_c$ and those derived from the fundamental parameters adopted in the literature for the spectral types studied in the present work (cf. Collins et al. 1991; Lang 1992; Townsend et al. 2004; Frémat et al. 2005; Huang & Gies 2006). Such a comparison is attempted in Fig. 1, where the $V_c$ parameters were calculated for masses and radii corresponding to the average MS spectral type-luminosity classes according to the different calibrations. Filled symbols and the respective $1\sigma_{V_c}$ dispersion bars represent the $V_c$ obtained from fundamental parameters given in the literature for the average MK spectral types (see Appendix A.2), while the open symbols correspond to our estimates of $V_c$. Figure 1 shows that, depending on the spectral type and luminosity class, our $V_c$ estimates can differ at most by 10 to 20 km s$^{-1}$ with those given by other authors.

Since in the BCD spectral classification system a given MK spectral type-luminosity class covers a curvilinear quadrilateral in the $(\lambda_1, D)$ diagram, our estimates of $V_c$ shown in Fig. 1, and given in Table 2, correspond to the centroids of these quadrilaterals, which do not necessarily coincide with those implied by the average MK spectral types defined by the stellar samples studied by other authors. In this table, and in Fig. 1, the spectral types are represented with the BCD index $S_{70}$ given in dex.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Filled symbols: average critical velocities $V_c$ and the related dispersions against the MK spectral type for luminosity classes V, IV, and III. These $V_c$ were calculated using the effective temperatures and bolometric luminosities from calibrations carried out by various authors. Open symbols: critical velocities for the average MK spectral types and luminosity classes that conform with the method used in this work to obtain the $V_c$ for individual stars. The MK spectral types are represented with the BCD parameter $S_{70}$ dex.</figcaption></figure>

<table><caption>**Table 2.** Critical velocities $V_c$ derived in this work for average MK spectral types and luminosity classes.</caption><thead><tr><th rowspan="2">Sp.T.</th><th rowspan="2">S<sub>70</sub><br>dex</th><th colspan="3">V<sub>c</sub> (km s<sup>-1</sup>)</th></tr><tr><th>V</th><th>IV</th><th>III</th></tr></thead><tbody><tr><td>O7</td><td>0.039</td><td>561</td><td>556</td><td>556</td></tr><tr><td>O9</td><td>0.052</td><td>542</td><td>538</td><td>517</td></tr><tr><td>B0</td><td>0.075</td><td>544</td><td>512</td><td>477</td></tr><tr><td>B1</td><td>0.104</td><td>530</td><td>494</td><td>423</td></tr><tr><td>B2</td><td>0.143</td><td>486</td><td>450</td><td>404</td></tr><tr><td>B3</td><td>0.181</td><td>461</td><td>418</td><td>380</td></tr><tr><td>B4</td><td>0.209</td><td>437</td><td>397</td><td>374</td></tr><tr><td>B5</td><td>0.238</td><td>438</td><td>381</td><td>362</td></tr><tr><td>B6</td><td>0.266</td><td>403</td><td>369</td><td>348</td></tr><tr><td>B7</td><td>0.297</td><td>393</td><td>360</td><td>341</td></tr><tr><td>B8</td><td>0.330</td><td>389</td><td>359</td><td>333</td></tr><tr><td>B9</td><td>0.384</td><td>386</td><td>355</td><td>326</td></tr><tr><td>A0</td><td>0.445</td><td>382</td><td>346</td><td>324</td></tr><tr><td>A1</td><td>0.491</td><td>375</td><td>343</td><td>322</td></tr></tbody></table>

**Notes.** The MK spectral types are represented with the BCD index $S_{70}$ dex.

## 3. Distribution of “true” rotational velocities

In this section we determine the distribution of ratios $V\sin i/V_c$, where the $V\sin i$ are corrected only for measurement uncertainties. This distribution is used in the rest of this work as the reference distribution of the data relative to rotational velocities. In a second step, we transform the distribution of projected ratios $V\sin i/V_c$ into a distribution of ratios $V/V_c$ of true rotational velocities. For this transformation we discuss two methods in
<figure><img src="image_3.png" /><figcaption>Fig. 2. a) Histogram giving the number N(Vsin i/Vc) of observed ratios v = Vsin i/Vc per class Δv studied in this work; b) smoothed density distribution Ψ(u) of projected velocity ratios u = Vsin i/Vc, where measurement uncertainties on u = Vsin i were taken into account. The error bars indicate the statistical uncertainty affecting the smoothed Ψ(u) distribution.</figcaption></figure>

which Stoeckley's corrections (Stoeckley 1968) of the observed Vsin i parameters are considered in different ways.

Since we shall deal with distributions corrected successively for a series of different effects, we adopt the notations v and Ψ(v) for velocity ratios and the corresponding distribution function as they are “before correction for some geometrical or physical effect”. The notations u and Φ(u) are used for velocity ratios and their distribution “after correction”, respectively. These notations can concern either the projected velocity ratios Vsin i/Vc or the true ratios V/Vc. We may then speak of v = Vsin i/Vc and Ψ(v) before correction for measurement uncertainties of the v values, which after correction become u = Vsin i/Vc and Φ(u), respectively. In the same way we have v = V/Vc and Ψ(v) before correcting the v ratios, for example for binary effects, which convert into u = V/Vc and Φ(u) once the correction is applied.

### 3.1. Distribution of rotational velocities corrected for measurement uncertainties

To obtain a first glance on the aspect of the distribution of apparent ratios v = Vsin i/Vc, we show in Fig. 2 the histogram constructed using the whole stellar sample regardless of the evolutionary state of the individual stars. The effect of the stellar evolution is minimized by considering Vsin i/Vc velocity ratios, where Vc is calculated consistently with the mass and evolutionary state of each object. Nonetheless, we discuss the possible evolutionary effects on the V/Vc ratios in Sect. 3.6. The class-steps of the histogram were established according to the bin-width optimization method by Shimazaki & Shinomoto (2007). We then obtained a smoothed version of the frequency density distribution of ratios Vsin i/Vc corrected for measurement uncertainties, which represents the reference distribution of the observed projected rotational velocities for the present work. This

<figure><img src="image_4.png" /><figcaption>Fig. 3. Blue line: relation between the observed (Vsin i)corr (ordinates) against the Vsin i corrected from Stoeckley's underestimation (abscissas). This relation was calculated for model He I 4471 line in stars with pnrc parameters T<sub>eff</sub> = 21 000 K and log g = 4.0, inclination angle i = 70°, and η ratios ranging from 0 to 1.0. The red bars indicate Stoeckley's correction Σ for an arbitrary near critical (Vsin i)corr parameter. The η values indicated in the figure are for the actual (Vsin i)<sub>Σ</sub> = V(η)sin i parameters in abscissas.</figcaption></figure>

distribution was established using kernel estimators (Bowman & Azzalini 1997), where each observed parameter v = Vsin i/Vc was represented by a Gaussian distribution, whose dispersion is given by the standard deviation of individual estimates. The distribution Ψ(u) thus obtained of ratios u = (Vsin i/Vc)corr, i.e., observed ratios v corrected for measurement uncertainties, is represented with a blue curve in Fig. 2. The error bars indicate the statistical uncertainty affecting the determination of the so-called smoothed distribution Ψ(u).

### 3.2. The Stoeckley effect

In rapid rotators, the emitted bolometric flux weakens from the pole to the equator as a function of the surface effective gravity (von Zeipel 1924; Lucy 1967; Espinosa Lara & Rieutord 2011). The contribution of the radiation to the total λ-dependent flux in a spectral line broadened by the rotation is thus less effective from the equatorial regions, which consequently translates into an underestimation of the Vsin i parameter. As this effect was first studied by Stoeckley (1968), in what follows we refer to it as the Stoeckley effect. The expected value of (Vsin i)<sub>Σ</sub> corrected from the Stoeckley effect can be written as

$$ (V\sin i)_\Sigma = (V\sin i)_{\text{corr}}^{\text{obs}} + \Sigma(\eta, i, M, t), \quad (7) $$

where (Vsin i)corr is the observed projected rotational velocity corrected for measurement uncertainties; Σ is the correction for the rotationally induced underestimation that hereafter, and depending on the circumstances, we call Stoeckley's correction or Stoeckley's underestimation. A schematic representation of the relation between the true Vsin i and the observed one is shown in Fig. 3.

After Stoeckley (1968), a number of authors have calculated Σ, in particular Townsend et al. (2004), Cranmer (2005), and Frémat et al. (2005), who calculated this correction strictly in
the frame of rigid rotation with the classical formulation of the GD effect by von Zeipel (1924). Frémat et al. (2005) calculated spectral lines in the visible spectral range in great detail using non-LTE model atmospheres of massive and intermediate-mass stars, mainly for He I 4471 and Mg II 4481 lines.

It is worth mentioning, however, that in spite of the colatitude $\theta$-dependent GD exponent $\beta_1$ (Espinosa Lara & Rieutord 2011; Rieutord 2016; Zorec et al. 2016), the effective temperature is still written as (von Zeipel 1924; Tassoul 1978)

$$T_{\text{eff}}^4(\theta) = \text{const.} g_{\text{eff}}^{\beta_1}(\theta), \quad (8)$$

with $\beta_1 = 1$. Other formulations of this relation for conservative and nonconservative rotation laws, produce $\beta_1 < 1$ (e.g., von Zeipel 1924; Lucy 1967; Smith & Worley 1974; Kippenhahn 1977; Hadrava 1992; Maeder 1999, 2009; Lovekin et al. 2006; Gillich et al. 2008; Dall & Sbordone 2011; Espinosa Lara & Rieutord 2011; Claret 2012). Moreover, because of the simple fact that $\beta_1$ is colatitude $\theta$-dependent, an observed $\beta_1$ is necessarily a function of the aspect angle (Domiciano de Souza et al. 2014; Rieutord 2016; Zorec et al. 2016). Holding forcibly $\beta_1 = \text{const.}$ over the stellar surface, form Eqs. (4), (27), and (28) in Espinosa Lara & Rieutord (2011) and the condition that for rigid rotation the ratio between the equatorial and polar radii is $R_e/R_p = 1 + \eta/2$, readily yields

$$\beta_1 = 1 - \frac{2}{3} \left[ \frac{\ln(1 - \eta) + \eta(1 + \eta/2)^{-3}}{\ln(1 - \eta) - 2 \ln(1 - \eta/2)} \right], \quad (9)$$

which shows that this exponent varies from $\beta_1 = 1$ at $\eta = 0$ to $\beta_1 = 1/3$ when $\eta = 1$. The function $\beta = \beta_1/4$ against the stellar flattening $\epsilon = 1 - R_p/R_e = \eta/(2 + \eta)$ is shown in Espinosa Lara & Rieutord (2011) and Domiciano de Souza et al. (2014). According to this dependence of $\beta_1$ with $\eta$, the contrast between the emitted polar and equatorial radiative fluxes reduces significantly when $\eta \to 1$ as compared to predictions obtained with $\beta_1 = 1$. We adopted $\beta_1 = 1$ on purpose, however, to produce the highest possible contrasts of effective temperatures in the stellar surface and maximize in this way the correction $\Sigma(\eta, i, M, t)$. This favors artificially inferring the highest possible values of $V$ in the $V\sin i$ parameter. Following these assumptions, we recalculated $\Sigma$ with FASTROT (Frémat et al. 2005) imposing rigid rotation and, for angular velocity, ratios ranging from $\Omega/\Omega_c = 0.0$ to $\Omega/\Omega_c = 0.9999$ ($0 \le \eta \le 0.985$).

### 3.3. Deconvolution methods

When the inclination angle of the stellar rotational axis is assumed to be distributed at random, the probability of occurrence of the inclination angle between $i$ and $i + di$ is given by $P(i)di = \sin i di$. To derive the distribution $\Phi(u)$ of ratios of true velocities $u = V/V_c$ from the distribution $\Psi(v)$ of ratios of apparent velocities $v = V\sin i/V_c$, we can adapt the rules given in Appendix D by writing $v = u \sin i$. The variables $(x, y, z)$ in Appendix D then translate into $(u, i, v)$ and the functions $Z(z)$ and $X(x)$ into $\Psi(v)$ and $\Phi(u)$, respectively. Accordingly, $dz/(\partial z/\partial y)$ becomes $di = dv/(u \cos i)$ and Eq. (D.3) transforms into the known Abel-like integral equation, first discussed by Kuiper (1935) and Chandrasekhar & Münch (1950)

$$\Psi(v) = v \int_{v}^{\infty} \frac{\Phi(u)}{u(u^2 - v^2)^{1/2}} du. \quad (10)$$

Cranmer (2005) notes that when we have to deal with the corrected velocity ratio $v = u \sin i + \Sigma(\eta, i)$, Eq. (10) does not take

into consideration the effect of the term $\Sigma$, and the correspondence between $\Phi(u)$ and $\Psi(v)$ is no longer unique.

So, if we take this observation by Cranmer into account, we guess that at least two different ways can be envisioned to derive $\Phi(u)$ from $\Psi(v)$: a) using the method introduced by Cranmer (2005), where $\Phi(u)$ is given a generic analytic form and each ratio $v$ is corrected for the corresponding GD effect according to all possible values of $u$ and $i$ to form the same $v$; b) using the Lucy-Richardson method (Lucy 1974; Richardson 1972) to solve Eq. (10), where $v$ of each star is corrected in advance for Stoeckley's underestimation.

In Sects. 3.4 and 3.5 we apply both methods in turn and discuss the pros and cons of their use.

### 3.4. Cranmer method

In the method originally introduced by Cranmer (2005), the function $\Phi(u)^2$ is represented by a truncated linear function (trapezoid),

$$\bar{\Phi}(u) = p \times u + q; \quad u_m \le u \le u_M, \quad (11)$$

where four parameters have to be determined: the constants $p$ and $q$, and the two extreme values [$u_m, u_M$] over which is defined $\bar{\Phi}(u)$, so that $\bar{\Phi}(u) \neq 0$ when $u_m \le u \le u_M$, and $\bar{\Phi}(u) = 0$ outside the interval. Actually, only three of them are iterated, since $q$ is determined automatically through the normalization condition $\int_{u_m}^{u_M} \bar{\Phi}(u) du = 1$. These parameters are inferred by fitting the function $\Psi(v)$ with a series of Monte Carlo trials using $N = 10^5$ model stars. In this method, the basic requirement of Cranmer is solved by assigning Stoeckley's correction to each built projected velocity ratio $v = u \sin i$, where the quantity $u$ belongs to a defined basis function and $i$ is assigned at random. The ratios $v$ are treated with an imposed uncertainty.

We tested the Cranmer method, but introduce small formal changes as follows:

(1) The simulations were tested using the distribution of projected rotational velocities corrected in advance for measurement uncertainties, for example, the distribution obtained in Sect. 3.1. According to this, the ratios $v = u \sin i$ simulated in the Monte Carlo trials are considered free from measurement uncertainties.

(2) In Appendix F it is shown that within some limiting conditions we can choose many different analytical expressions for $\Phi(u)$, which all can be considered statistically equivalent. We then replaced the above trapezoidal form by the following function,

$$\bar{\Phi}(u) = Au^a \exp[-b(u - c)^2], \quad (12)$$

which also depends on four free parameters: $A$, $a$, $b$, and $c$, where $A$ is fixed by the normalization condition $\int_{0}^{1} \Phi(u) du = 1$. This analytical form was chosen because it resembles the $\Phi(u)$ function determined by the Lucy-Richardson deconvolution method of data analyzed in the present work (see Sect. 3.5).

(3) Stoeckley's correction were assigned to the $v$ parameters not only as a function of the $u$- and $i$-Monte Carlo simulated parameters, but also depending on ($T_{\text{eff}}$, log $g$). The effective temperatures and gravities used are the average $T_{\text{eff}}$ and log $g$.

<footnote n="2">² The notation $\Phi(u)$ is used for the actual distribution, while $\bar{\Phi}(u)$ indicates its analytical/auxiliary representation.</footnote>
<table><caption>**Table 3.** Solution parameters (a, b, c, A) estimated using the Cranmer method for the function $\Phi(u)$ defined in Eq. (12).</caption><thead><tr><th></th><th>a</th><th>b</th><th>c</th><th>A</th><th>u<sub>mode</sub></th></tr></thead><tbody><tr><td>$\bar{\Phi}(u)_{Cr1}$</td><td>0.000</td><td>44.179</td><td>0.690</td><td>3.722</td><td>0.6900</td></tr><tr><td>$\bar{\Phi}(u)_{Cr2}$</td><td>8.809</td><td>23.517</td><td>0.420</td><td>475.788</td><td>0.6912</td></tr></tbody></table>

values calculated from the apparent parameters of the studied stellar sample distributed in 16 equidistant groups. For each tested triplet (a, b, c), the fraction of model stars with the same ($T_{\text{eff}}$, log g) is determined by the fraction of real stars in the corresponding ($T_{\text{eff}}$, log g) group.

(4) The method was applied twice: (i) using a predefined basis u-function; (ii) with a Monte Carlo sorting of the u parameters, using the probability function $\Phi(u)$ defined by the ongoing tested a, b, and c parameters. We noted that both methods lead to the same result.

(5) The N simulated $v=u \sin i$ parameters were gathered into an histogram whose class width, determined as in Sect. 3.1, was adopted as the dispersion of the Gaussian kernels needed to obtain a smoothed version of the predicted frequency density distribution of apparent velocity ratios. As in Cranmer's method, the quality of the fit of the predicted frequency with the reference distribution $\Psi(v)$ was controlled with the $\chi^2$ test.

Curiously, two sets of parameters (a, b, c) for $\bar{\Phi}(u)$ given by Eq. (12) were obtained that produce the same 'best' fit of $\Psi(v)$ according to the $\chi^2$ test. We call these solutions $\bar{\Psi}(u)_{Cr1}$ and $\bar{\Phi}(u)_{Cr1}$, respectively. They are shown in Fig. 4 and the parameters characterizing them are given in Table 3. For comparison sake, Fig. 4 also shows the reference distribution $\Psi(v)$ obtained in Sect. 3.1 (see Fig. 2 curve b). Nevertheless, both solutions $\bar{\Phi}(u)_{Cr1}$ and $\bar{\Phi}(u)_{Cr2}$ are probably equivalent within the uncertainties that affect $\Psi(v)$ shown in Fig. 2 (curve b).

### 3.5. Lucy-Richardson deconvolution method

In the original Lucy-Richardson deconvolution method, the solution for $\Phi(u)$ is obtained from Eq. (10) with a subsidiary equation constructed using a Bayesian equivalence principle. According to this principle the probability of occurrence of u and v can be written either as $\Phi(u)P(u|v)$ or $\Psi(v)P(v|u)$, where for a value v the probability of the occurrence of u is represented with $P(u|v)$, presently given by the factor multiplying $\Phi(u)$ in the integrand of Eq. (12). This equivalence translates into

$$ \Phi^{r+1}(u) = \Phi^r(u) \int_{\mathcal{V}(u)} \left[ \frac{\Psi(v)}{\Psi^r(v)} \right] P(v|u) dv, \quad (13) $$

which together with Eq. (10) enables iteration of $\Phi(u)$ and $\Psi(v)$. In Eq. (13), $\Phi^r(u)$ and $\Psi^r(v)$ are the r-iterated estimates of the respective distributions, and $\mathcal{V}(u)$ represents the integration interval dependent on u. The optimal number of required iterations was controlled through the Kolmogorov-Smirnov test over $\Phi(u)$ and on the recalculated $\Psi(v)$ compared with the original or input apparent distribution.

To account for the main Cranmer demand, i.e., to consider Stoeckley's correction consistently, the application of the Lucy-Richardson method becomes possible when each apparent ratio v is corrected for Stoeckley's underestimation in advance of the iteration process. In the next section we detail our procedure

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Diagram showing the two functions $\bar{\Phi}(u)$ that lead to the same $\chi^2$; $\bar{\Phi}(u = V/V_c)_{Cr1}$ (red full line) and $\bar{\Phi}(u = V/V_c)_{Cr2}$ (red dashed line), obtained with the Cranmer method and the analytical form given in Eq. (12). The reference distribution $\Psi(v)$ of ratios $v = V \sin i / V_c$ corrected for observational uncertainties is also shown (blue dashed curve).</figcaption></figure>

to estimate $\Sigma(\eta, i, M, t)$ and derive $v = V \sin i / V_c$ corrected for Stoeckley's underestimation.

#### 3.5.1. Correction $\Sigma(\eta, i, M, t)$ as a function of pnrc parameters

To estimate $\Sigma(\eta, i, M, t)$ we assume that stars are rigid rotators and use Eq. (8) with $\beta_1 = 1$. Stoeckley's correction is then a function of the stellar mass $M/M_\odot$, stellar age $t/t_{\text{MS}}$ ($t_{\text{MS}}$ is the time a rotating star spends in the main-sequence phase), the ratio of the centrifugal to gravitational acceleration at the equator $\eta$ (see Eq. (1)), and inclination $i$ of the rotation axis. To determine these quantities, the following form is given to the stellar apparent fundamental parameters:

$$ \begin{align} T_{\text{eff}}^{\text{app}} &= T_{\text{eff}}^{\text{pnrc}}(M, t) C_{\text{T}}(M, t, \eta, i) \nonumber \\ g_{\text{eff}}^{\text{app}} &= g_{\text{eff}}^{\text{pnrc}}(M, t) C_{\text{G}}(M, t, \eta, i) \nonumber \\ L_{\text{app}}^{\text{app}} &= L_{\text{pnrc}}^{\text{app}}(M, t) C_{\text{L}}(M, t, \eta, i) \tag{14} \end{align} $$

$$ \frac{(V \sin i)_{\text{app}}}{V_c(M, t)} = \left[ \frac{\eta}{R_e(M, t, \eta)/R_c(M, t)} \right]^{1/2} \sin i - \frac{\Sigma(M, t, \eta, i)}{V_c(M, t)}, $$

where $R_e(M, t, \eta)$ and $R_c(M, t)$ are the actual and critical stellar equatorial radii, which are determined using our 2D models of rigidly rotating stars (Zorec et al. 2011; Zorec & Royer 2012). The left side of Eq. (3) is associated with the observed (apparent) fundamental parameters determined in Sect. 2.2. On the right side of Eq. (3), $T_{\text{eff}}^{\text{pnrc}}(M, t)$, $g_{\text{eff}}^{\text{pnrc}}(M, t)$, $L_{\text{pnrc}}^{\text{app}}(M, t)$ are the pnrc effective temperature, surface gravity and bolometric luminosity introduced in Sect. 2.2. The functions $C_{\text{T}}(M, t, \eta, i)$, $C_{\text{G}}(M, t, \eta, i)$, and $C_{\text{L}}(M, t, \eta, i)$ carry all the information relative to the geometrical deformation of the rotating star and of its GD over the observed hemisphere.

To determine $C_{\text{T}}(M, t, \eta, i)$ and $C_{\text{G}}(M, t, \eta, i)$, we calculated synthetic spectra in the $\lambda\lambda 4000-4500$ Å wavelength interval as
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** a) Function $C_T(M, t, \eta, i)$; b) function $C_L(M, t, \eta, i)$; c) function $C_G(M, t, \eta, i)$. The pnrc parameters used here are $T_{\text{eff}}^{\text{pnrc}} = 20\,000$, $\log g_{\text{pnrc}} = 3.8$ and inclination angles in the interval $0 \le i \le \pi/2$ at steps $\Delta i = 10^\circ$.</figcaption></figure>

emitted by our 2D models of rigid rotators seen at several aspect angles $i$. The effective temperature and surface gravity of the rotationless counterparts of these models define the respective $T_{\text{eff}}^{\text{pnrc}}(M, t)$ and $g_{\text{eff}}^{\text{pnrc}}(M, t)$ parameters. The effective temperatures and surface gravities of classic nonrotating, plane-parallel model atmospheres, which fit the above synthetic spectra of rotating stars, define the theoretical $T_{\text{eff}}^{\text{app}}$ and $g_{\text{eff}}^{\text{app}}$ parameters; this enables us to calculate the functions $C_T(M, t, \eta, i) = T_{\text{eff}}^{\text{app}}/T_{\text{eff}}^{\text{pnrc}}$ and $C_G(M, t, \eta, i) = g_{\text{eff}}^{\text{app}}/g_{\text{eff}}^{\text{pnrc}}$. The fit of spectra was performed with the CERN-MINUIT minimization package.

The function $C_L(M, t, \eta, i)$ was determined employing a similar formulation as in Georgy et al. (2014), except that in Eq. (8) we maintained $\beta_1 = 1$ on purpose for all values of $\eta$ to produce the largest possible Stoeckley’s corrections $\Sigma(M, t, \eta, i)$. Figure 5 shows the functions $C_T(M, t, \eta, i)$, $C_L(M, t, \eta, i)$ and $C_G(M, t, \eta, i)$ calculated for an object with mass $M/M_\odot = 8.2$ and age $t/t_{\text{MS}} = 0.9$, i.e., pnrc parameters $T_{\text{eff}}^{\text{pnrc}} = 20\,000$ K and $\log g^{\text{pnrc}} = 3.8$.

As a result of the known “mass-lowering effect” induced by the rotation (Milne 1923; Sackmann 1970; Bodenheimer 1971; Tuominen 1972; Clement 1979), according to which a rotating star behaves as another object with lower mass and, thus, its bolometric luminosity is lower the higher its rotation, these functions depend sensitively on the stellar mass and age. Hence, the coefficients $C_T$, $C_L$, and $C_G$ calculated for the present work do not exactly reproduce the behavior of the absolute magnitude $M_V$ against the color ($B-V$) as a function of $\eta$ and the inclination $i$ obtained by Collins et al. (1991). In fact, the magnitude $M_V$ determined by Collins et al. (1991) remains brighter than that of a nonrotating counterpart even at $i \gtrsim 60^\circ$ as $\eta \to 1$ for spectral types hotter than B3. This is because of two important simplifications introduced in their calculations: (1) the “mass-lowering effect” was neglected; and (2) the polar radius $R_p$ is maintained unchanged, even though its variation with $\eta$ introduces non-negligible effects on the magnitudes as $i \gtrsim 50^\circ$. The changes of the bolometric luminosity used in the present work are detailed in Frémat et al. (2005), while the variation of the $R_e$ and $R_p$ radii are given in Zorec et al. (2011) and Zorec & Royer (2012). We recall that for $M/M_\odot = 3$ the lowering in the bolometric magnitude ranges from $\Delta M_{\text{bol}} = +0.12$ to +0.22 mag as $\eta$ changes from 0.3 to 1.0, while for $M/M_\odot = 20$ we have $\Delta M_{\text{bol}} = +0.08$ mag at $\eta = 0.3$ and $\Delta M_{\text{bol}} = +0.12$ mag at $\eta = 1.0$. The respective changes $\Delta M_V$ are slightly smaller because of the additional bolometric corrections. In any case, the aspect-angle dependence of $M_V$ with $i \gtrsim 50^\circ$ and $\eta$ are much larger than those predicted by Collins et al. (1991).

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Distribution function $\Phi(u)_{\text{LR}}$ obtained with the Lucy-Richardson method (blue curve) with the corresponding statistical uncertainties. Both distributions $\bar{\Phi}(u)_{\text{Cr1,2}}$ of true rotational velocities derived using the Cranmer method (red curves) are also superimposed.</figcaption></figure>

The solution of Eq. (3) automatically determines the $\Sigma(M, t, \eta, i)$ that corrects the apparent $V\sin i$ of each star for the GD effect in a consistent way with the stellar fundamental parameters $M/M_\odot$, $t/t_{\text{MS}}$, $\eta$ and $i$. Details on the method used to solve the system of equations in Eq. (3) are given in Sect. E. The pnrc effective temperatures, gravities, bolometric luminosities, ratios $\eta$, and inclinations $i$ obtained for the program Be stars are listed in Table 4. We did not consider the stars HD 49131, HD 50737, HD 112091, HD 166566, HD 171054, and HD 330950 in this work because we obtained uncertain pnrc fundamental parameters, which situate then far below the ZAMS.

We then determined the smoothed distribution of projected velocity ratios $v = [V\sin i + \Sigma(M, t, \eta, i)]/V_c$ corrected for measurement uncertainties of $V\sin i$, using the same procedure as in Sect. 3.1 that by iteration of Eqs. (10) and (13) enabled us to derive the distribution $\Phi(u)_{\text{LR}}$ of ratios $u = V/V_c$, which is shown in
<table><caption>**Table 4.** *pnrc* fundamental parameters of the studied Be stars.</caption><thead><tr><th rowspan="2">HD</th><th rowspan="2">T<sup>pnrc</sup><sub>eff</sub></th><th rowspan="2">± σ<sub>T<sub>eff</sub></sub></th><th>log g<sub>pnrc</sub></th><th>log L<sub>pnrc</sub>/L<sub>☉</sub></th><th>(V sin i)<sub>pnrc</sub></th><th>V<sub>c</sub> ± σ<sub>V<sub>c</sub></sub></th><th>M/M<sub>☉</sub> ± σ<sub>M</sub></th><th>η ± σ<sub>η</sub></th><th>i ± σ<sub>i</sub></th></tr><tr><th>± σ<sub>log g</sub></th><th>± σ<sub>log L</sub></th><th>± σ<sub>Vsin i</sub></th><th>km s<sup>-1</sup></th><th>km s<sup>-1</sup></th><th>deg</th></tr><tr><th></th><th>K</th><th></th><th>dex</th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>144</td><td>13 300 ± 180</td><td>3.81 ± 0.10</td><td>2.912 ± 0.028</td><td>140 ± 11</td><td>318 ± 15</td><td>4.59 ± 0.23</td><td>0.55 ± 0.08</td><td>38 ± 9</td></tr><tr><td>4180</td><td>17 400 ± 500</td><td>3.70 ± 0.31</td><td>3.514 ± 0.057</td><td>215 ± 18</td><td>357 ± 19</td><td>6.59 ± 0.33</td><td>0.47 ± 0.14</td><td>58 ± 14</td></tr><tr><td>5394</td><td>40 200 ± 2320</td><td>4.25 ± 0.34</td><td>5.274 ± 0.084</td><td>441 ± 51</td><td>677 ± 59</td><td>31.69 ± 2.01</td><td>0.52 ± 0.13</td><td>65 ± 16</td></tr><tr><td>330 950</td><td>30 440 ± 3650</td><td>4.54 ± 0.39</td><td>2.475 ± 0.121</td><td>60 ± 6</td><td>531 ± 91</td><td>4.94 ± 0.55</td><td>0.47 ± 0.11</td><td>8 ± 5</td></tr></tbody></table>

**Notes.** This table is given in full at the CDS.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Distributions $\Phi(u)_{LR}$ obtained with Stoeckley's corrections $\Sigma(\eta, \beta_1 = 1)$ calculated for different values of $\eta$: $\eta = \eta_*$ (derived in Sect. 3.5.1); $\eta = \eta_{min}$; $\eta = 1.0$ and $\eta = 0.0$. For comparison the smoothed distribution of ratios of *pnrc* velocities $V\sin i/V_c\sin i$ (crosses) is also shown.</figcaption></figure>

Fig. 6 (blue curve). This figure shows the statistical uncertainties of $\Phi(u)_{LR}$ and, for comparison, the two distributions $\Phi(u)_{Cr1,2}$ derived with Cranmer's method (red curves).

Since $\Phi(u)_{LR}$ depends on Stoeckley's corrections $\Sigma$ obtained by solving Eq. (3), it is worth asking what uncertainties affect its determination that rely on the estimate of $\eta$. To this end we have tested three extreme cases: (a) all stars have the lower possible rotational parameter $\eta_{min}$, i.e., the value derived when we consider that $V(\eta_{min}) = V\sin i$; (b) Stoeckley's corrections estimated for $\eta = 1$, i.e., $V(\eta = 1) = V_c$; (c) neglecting Stoeckley's correction, i.e., $\Sigma(M, t, \eta = 0, i) = 0$ in all stars. The distribution functions $\Phi(u)_{LR}$ thus obtained are shown in Fig. 7. Also, the solution that comes from $\eta = \eta_*$ derived with Eq. (3) (blue full line) is also shown in this section. As expected, approximation (a) produces a slight excess of low rotators (dashed red line); solution (b) is characterized by an excess of rapid rotators (red full line); (c) the distribution is shifted to smaller values of $u = V/V_c$ (blue long-dashed line).

Another test of consistency of both the *pnrc* parameters used to infer Stoeckley's correction and of $\Phi(u)_{LR}$ can be performed by calculating just the smoothed distribution $\Phi(u_{pnrc})$ of the true

velocity ratios defined as $u_{pnrc} = [(V\sin i)_{pnrc}/V_c]/\sin i_{pnrc}$, and thus obtained independently of the Lucy-Richardson method. The smoothed distribution is established using the kernel-estimator method (Bowman & Azzalini 1997) employed in Sect. 3.1. Figure 7 shows this function as a black cross-marked curve. This function, despite the errors that may affect the solutions issued from Eq. (3), in particular the inclination angles, adjusts the distribution $\Phi(u)_{LR}$, which is obtained in Sect. 3.5 with Eqs. (10) and (13), and is widely within the statistical uncertainty limits.

As noted above, we used Eq. (8) with $\beta_1 = 1$ to simulate the GD effect with the purpose of having the largest possible mode of the $\Phi(V/V_c)$ distribution. Nonetheless, because in real stars not only $\beta_1 < 1$, but also $\beta_1 = \beta_1(\theta)$ (Espinosa Lara & Rieutord 2011; Zorec et al. 2016), it is expected that the actual mode of $\Phi(V/V_c)$ be displaced to a even lower value that found here, i.e., $(V/V_c)_{mode} < 0.7$.

If the inclination angles were not distributed at random, $\Phi(u)_{LR}$ could suffer from an additional skewness, which is a function of that particular distribution of inclination angles. This issue is briefly discussed in Sect. 3.7.

### 3.5.2. Comments on the obtained distributions $\Phi(u)$ and adopted deconvolution method

Having obtained the distribution $\Phi(u)$ using two different methods, where both enforce the basic requirements for Stoeckley's correction of the apparent ratio $v = V\sin i/V_c$, we would like to know whether they can be considered equivalent to each other.

We can notice that the Lucy-Richardson method does not impose any restriction to the functional aspect of $\Phi(u)$, although it requires that Stoeckley's correction be applied before the iteration procedure, which is not a straightforward task.

The Cranmer method imposes the use of an analytical expression $\bar{\Phi}(u)$ meant to represent the actual function $\Phi(u)$, which necessarily depends on a limited number of free parameters. In this way, only a limited number of moments of the actual distribution $\Phi(u)$ can be accounted for. Moreover, when a deconvolution has to be carried out with kernels determined by complicated probability distributions, doubts may arise as to whether the chosen analytical function $\bar{\Phi}(u)$ may correctly account for the physical phenomena that shape the actual $\Phi(u)$ distribution.

Nonetheless, the application of the Cranmer method is simple and the Monte Carlo simulation can be significantly shortened when an approximation of the free parameters defining $\bar{\Phi}(u)$ is known in advance. In Appendix F we develop a simple analytical way to calculate these quantities for trapezoidal and triangular distributions $\bar{\Phi}(u)$. If the apparent variable $v$ is corrected in advance for Stoeckley's underestimation, the
<table><caption>**Table 5.** Parameters defining trapezoidal and triangular $\bar{\Phi}(u)$ derived according to the Cranmer method.</caption><thead><tr><th></th><th>u<sub>m</sub></th><th>u<sub>M</sub></th><th>p</th><th>q</th></tr></thead><tbody><tr><td>Φ̄(u)<sub>trapezoid</sub></td><td>0.468</td><td>0.900</td><td>1.856</td><td>1.044</td></tr><tr><td>Φ̄(u)<sub>triangle</sub></td><td>0.381</td><td>0.987</td><td>0.721</td><td>3.295</td></tr></tbody></table>

Notes. $u \in [u_m, u_M]$; $\bar{\Phi} = p.u + q$; $\bar{\Phi}_{max} = \bar{\Phi}(u_{max})$.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** a) Distribution function $\Phi_{LR}(u)$ and the associated auxiliary trapezoidal and triangular distributions $\bar{\Phi}(u)$. b) Cumulative distributions $C[\Phi(u)] \equiv \int_0^1 \Phi(u)du$ of the density distributions shown in a).</figcaption></figure>

analytical method produces the final values of the sought free parameters.

We compared the Cranmer and Lucy-Richardson methods using quantitative procedures. To this end, we have applied Cranmer's algorithm in full to our stellar sample using Eq. (12), trapezoidal, and triangular functions to represent $\bar{\Phi}(u)$. The free parameters thus obtained are given in Table 5, and the respective functions are shown in Fig. 8a. We then applied two different tests: (i) comparison of the first four moments of distributions and (ii) the Kolmogorov-Smirnov test.

### (i) Moments of distributions

The first four moments are usually employed to control four basic characteristics of distributions: the average value of the studied random variable; the variance, or dispersion of the studied random variable; the skewness, or symmetry of distributions; and the kurtosis, or flattening of distributions. In Table 6 are given the moments $\langle u^n \rangle$, with $n = 1, ..., 4$, of functions $\Phi(u)_{LR}$, $\Phi(u)_{Cr1,2}$, $\Phi(u)_{trapezoid}$, and $\Phi(u)_{triangle}$, as well as their respective modes ($u_{max}$).

From Fig. 6 and the quantities given in Tables 5 and 6 we can then conclude:

(1) The distributions, $\Phi(u)_{Cr1,2}$ and $\Phi(u)_{LR}$ are defined over the same interval of $u$ ratios and their modes are fairly similar.

(2) The Cranmer solution $\Phi(u)_{Cr2}$ is closer to $\Phi(u)_{LR}$. The approach $\Phi(u)_{Cr1}$ suggests a relative lower number of rotators with $V/V_c \ge 0.8$ and a relative higher number of rotators, but it has nearly the same mode of velocities. This may be due to a compensation effect related to the use of an imposed analytical form for $\Phi(u)$.

<table><caption>**Table 6.** First four moments of the function $\Phi(u)$ derived according to different approaches used: Lucy-Richardson; Cranmer with pseudo-Maxwellian, trapezoidal and triangular functions.</caption><thead><tr><th>Methods</th><th>⟨u<sup>1</sup>⟩</th><th>⟨u<sup>2</sup>⟩</th><th>⟨u<sup>3</sup>⟩</th><th>⟨u<sup>4</sup>⟩</th><th>u<sub>max</sub></th><th>p%</th></tr></thead><tbody><tr><td>Φ(u)<sub>LR</sub></td><td>0.697</td><td>0.501</td><td>0.370</td><td>0.281</td><td>0.684</td><td>—</td></tr><tr><td>Φ̄(u)<sub>Cr1</sub></td><td>0.689</td><td>0.487</td><td>0.351</td><td>0.258</td><td>0.690</td><td>60.0</td></tr><tr><td>Φ̄(u)<sub>Cr2</sub></td><td>0.695</td><td>0.497</td><td>0.365</td><td>0.275</td><td>0.691</td><td>99.9</td></tr><tr><td>Φ̄(u)<sub>trapezoid</sub></td><td>0.696</td><td>0.500</td><td>0.370</td><td>0.279</td><td>0.735*</td><td>20.5</td></tr><tr><td>Φ̄(u)<sub>triangle</sub></td><td>0.696</td><td>0.500</td><td>0.370</td><td>0.279</td><td>0.721</td><td>99.4</td></tr></tbody></table>

Notes. The last column gives the confidential level at which the Cranmer distributions can be considered parent to the Lucy-Richardson distribution. (*) The parameter given here is the median of the distribution.

(3) The first three moments produced by the trapezoidal and triangular auxiliary distributions using the Cranmer algorithm in full are closer to those of $\Phi(u)_{LR}$ than those due to pseudo-Maxwellian auxiliary distributions.

(4) Different analytical forms of $\bar{\Phi}(u)$ are in principle able to reproduce distributions that reassemble the observed $\Phi(v)$, i.e., having nearly the same first four moments, even though the same probability $P(u|v)$ was employed in all cases.

We can then wonder whether imposed analytic functions $\bar{\Phi}(u)$ may reliably describe the properties of the actual and sought distribution $\Phi(u)$.

### (ii) The Kolmogorov-Smirnov test

To simplify the comparisons, we compared the distributions $\Phi(u)_{Cr1,2}$ and $\bar{\Phi}(u)_{trapez,triang}$ with $\Phi(u)_{LR}$. The cumulative distributions $C[\Phi(u)_{Cr1,2}]$, $C[\bar{\Phi}(u)_{trapezoid}]$ and $C[\bar{\Phi}(u)_{triangle}]$ are shown in Fig. 8b; the notation means $C(\Phi) = \int_0^1 \Phi(u) du$. The last column of Table 6 gives the probability at which the analytical functions $\bar{\Phi}(u)$ imposed here can be considered parent to $\Phi(u)_{LR}$.

According to the Kolmogorov-Smirnov test, we see that the Cranmer solutions $\Phi(u)_{Cr1,2}$ with pseudo-Maxwellian functions can be considered parent distributions of $\Phi(u)_{LR}$, respectively, within 60.0% and 99.9% levels of confidence, while $\bar{\Phi}(u)_{trapezoid}$ and $\bar{\Phi}(u)_{triangle}$ are at 20.5% and 99.4% confidence levels, respectively. Also, the use of trapezoidal distributions to discuss the velocity distributions of Be stars is less appropriate.

So, according to tests (i) and (ii), it seems sensible to use the Lucy-Richardson deconvolution method because it avoids imposing analytical behaviors for $\Phi(u)$. Moreover, the complexity of some corrections applied in the next sections, makes the use of analytical auxiliary functions to account for the inherent subtleties that may characterize $\Phi(u)$ rather uncertain and complicated.

## 3.6. Effects of stellar evolution on the distribution of $V/V_c$ ratios

In general, a stellar sample can be considered statistically reliable if it is formed by a rather large number of objects, which all have roughly the same masses and ages. When this is not so, however, the effects assumed by the differences in masses can be somewhat suppressed using normalized apparent rotational
<table><caption>**Table 7.** Main-sequence evolutionary subphases, number of program stars in each group, the corresponding average fractional times $\langle t/t_{\text{MS}} \rangle$ with their dispersions, modes $(V/V_c)_{\text{mode}}$ of the distributions, and the respective characteristic dispersions.</caption><thead><tr><th>Time interval</th><th>N<sub>*</sub> &lt; t/t<sub>MS</sub> &gt; &sigma;<sub>t/t<sub>MS</sub></sub> (V/V<sub>c</sub>)<sub>mode</sub> &plusmn; &sigma;<sub>V/V<sub>c</sub></sub></th><th></th><th></th></tr></thead><tbody><tr><td>1 0.00 &le; t/t<sub>MS</sub> &le; 0.61</td><td>77</td><td>0.25 &plusmn; 0.23</td><td>0.64 &plusmn; 0.12</td></tr><tr><td>2 0.61 &le; t/t<sub>MS</sub> &le; 0.90</td><td>79</td><td>0.77 &plusmn; 0.09</td><td>0.76 &plusmn; 0.10</td></tr><tr><td>3 t/t<sub>MS</sub> &ge; 0.90</td><td>77</td><td>1.01 &plusmn; 0.03</td><td>0.77 &plusmn; 0.15</td></tr><tr><td>4 0.00 &le; t/t<sub>MS</sub> &le; 0.33</td><td>44</td><td>0.06 &plusmn; 0.09</td><td>0.63 &plusmn; 0.09</td></tr><tr><td>5 0.33 &le; t/t<sub>MS</sub> &le; 0.66</td><td>48</td><td>0.54 &plusmn; 0.10</td><td>0.72 &plusmn; 0.15</td></tr><tr><td>6 t/t<sub>MS</sub> &ge; 0.66</td><td>141</td><td>0.91 &plusmn; 0.11</td><td>0.76 &plusmn; 0.13</td></tr><tr><td>7 0.66 &le; t/t<sub>MS</sub> &le; 1.00</td><td>85</td><td>0.84 &plusmn; 0.09</td><td>0.76 &plusmn; 0.10</td></tr></tbody></table>

velocities $V \sin i/V_c$. Nevertheless, the ratio $V/V_c$ depends on stellar evolutionary effects, which are controlled by the rotation itself (cf. Maeder 2009). It is then difficult to establish a statistically sample of Be stars where the effects of the evolution are entirely obliterated.

In this study we considered all Be stars forming a unique group, although they are in different MS evolutionary phases. It is however important to inquire which statistical tendencies can be brought out from rotation-velocity distributions of stars in different evolutionary time intervals, in particular on the evolution of the ratio $V/V_c$. For early-type stars, few studies have attempted to tackle this issue by putting stars into luminosity-class groups (Balona 1975; Zorec 1986; Zorec et al. 1988, 1990, 2004; Yudin 2001; Cranmer 2005). When dealing with rapid rotators, however, the MK luminosity classes are apparent stellar properties whose relation with the actual evolutionary characteristics is not straightforward. We then divided our Be-star sample into three groups of different fractional age $t/t_{\text{MS}}$-intervals: early MS period, intermediate MS period, and late MS period, using the $t/t_{\text{MS}}$ determined in Sect. 3.5.1. This division is marred by uncertainties. Adopting equal $t/t_{\text{MS}}$-intervals, the number of stars in each of these groups is very different; in particular, there are some which can be considered statistically ‘underpopulated’ groups. On the contrary, if the same number of stars is demanded for each group, the $t/t_{\text{MS}}$-intervals become very unequal. Some of them are quite large in duration, which cannot then single out characteristics proper to a given short evolutionary period. Table 7 gives the characteristics of some Be-star groups we could form. As a number of stars have $t/t_{\text{MS}} \ge 1$, we also established a group with objects in $0.66 \le t/t_{\text{MS}} \le 1.00$.

The distributions of $\Phi(u)$ velocity ratios $u = V/V_c$ obtained for these stellar groups are shown in Fig. 9, where the area under each curve is proportional to the number of stars in the group. It is thus clear that a significantly larger number of objects would be necessary to draw more precise results on the evolution of the ratio $V/V_c$. However, these distributions confirm the known fact that $V/V_c$ becomes globally larger as stars evolve from early MS evolutionary phases to later evolutionary phases, whether or not they are rapid rotators (Yudin 2001; Zorec et al. 2004; Cranmer 2005). Nevertheless, Fig. 9 uncovers another property that becomes evident when we compare the $V/V_c$-distributions of intermediate and late MS evolutionary periods, which is probably specific to Be stars: once the mode $(V/V_c)_{\text{mode}} \approx 0.77$ is attained at some intermediate evolutionary phase in the MS, it will stay unchanged up to the end of the MS evolutionary lifespan (see rows 2, 3, 6 and 7 in Table 7).

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** a) Distributions of velocity ratios $u = V/V_c$ of the program Be stars corresponding to three MS evolutionary subphases, which are indicated in the inlay. Each group has roughly the same number of stars. b) Distributions of velocity ratios $u = V/V_c$ corresponding to three MS evolutionary subphases of the same duration, which are indicated in the inlay. In each block, the normalization is made by adding the area of all three distributions.</figcaption></figure>

Studies of other effects that change the $V/V_c$-distributions of Be stars would certainly benefit from sample divisions by shorter evolutionary intervals than those that we outlined here. Nevertheless, on account of the statistical uncertainties just mentioned above, in the remainder of this work and in Paper II of this series (Zorec et al. 2016), we continue to treat the entire stellar sample forming a unique, single group.

The discussion of the variation of the ratio $V/V_c$ with the age can be extended a little by studying its dependence with the stellar mass. We distributed our stellar sample in age-mass groups with the goal of having the same average masses in the different age intervals. Unfortunately, this attempt ends up with groups where the number of objects are not the same and some of them are rather under-represented to the detriment of the quality of statistics. These groups and the corresponding average values of the *pncr* velocity ratios $V \sin i/V_c$ and $V/V_c$ are given in Table 8.

From Table 8 we can conclude that in each studied age interval, the dependence of the average of true velocity ratios with the stellar mass is marginal, if there is any. There is, however, a dependence with age in the sense that $V/V_c$ is higher as the age ratio $t/t_{\text{MS}} \to 1$ in all mass groups. Thus, this result does not confirm the claim by Cranmer (2005), which states that the lower limit of $V/V_c$ to form circumstellar disks in the hottest (more massive) Be stars is “well below 1”, while the lower limit approaches 1 in the less massive Be stars.

### 3.7. Effects on the distribution of ratios $V/V_c$ due to nonrandomly distributed inclinations

In Sects. 3.3–3.5 we assumed that the inclination angles are distributed at random, however, the studied sample may not agree with this hypothesis. To complete the discussion on the characteristics of distributions of true rotational velocities, we simulate cases where the inclination angle $i$ of the stellar rotational axis is not distributed at random and observe the induced effects on the distributions of ratios of true rotational velocities. Because there is not a unique way to formulate the lack of randomness, we assume the density probability of the occurrence of an inclination between $i$ and $i + di$ is given by

$$P_k(i) di = A_k \sin^{1+k} i \, di, \qquad (15)$$

where $k$ is a free parameter, while $A_k$ is the normalization constant. With $k > 0$, there are simulated frequencies favoring high
<table><caption>**Table 8.** Average *pncr* velocity ratios $\langle V \sin i / V_c \rangle$ and $\langle V / V_c \rangle$ for mass-agerage groups of the studied Be star sample.</caption><thead><tr><th rowspan="2"></th><th colspan="3">t/t<sub>MS</sub> ≤ 0.4</th><th colspan="3">0.4 &lt; t/t<sub>MS</sub> ≤ 0.8</th><th colspan="3">t/t<sub>MS</sub> &gt; 0.8</th></tr><tr><th>⟨M/M<sub>☉</sub>⟩</th><th>⟨t/t<sub>MS</sub>⟩</th><th>n</th><th>⟨M/M<sub>☉</sub>⟩</th><th>⟨t/t<sub>MS</sub>⟩</th><th>n</th><th>⟨M/M<sub>☉</sub>⟩</th><th>⟨t/t<sub>MS</sub>⟩</th><th>n</th></tr></thead><tbody><tr><td>M/M<sub>☉</sub> ≤ 7.5</td><td>5.2</td><td>0.09</td><td>20</td><td>5.0</td><td>0.62</td><td>29</td><td>5.6</td><td>0.97</td><td>34</td></tr><tr><td>7.5 &lt; M/M<sub>☉</sub> ≤ 15</td><td>10.1</td><td>0.11</td><td>20</td><td>10.7</td><td>0.63</td><td>40</td><td>10.5</td><td>0.96</td><td>40</td></tr><tr><td>15 &lt; M/M<sub>☉</sub> ≤ 25</td><td>19.0</td><td>0.12</td><td>7</td><td>18.4</td><td>0.75</td><td>9</td><td>18.2</td><td>0.98</td><td>15</td></tr><tr><td>M/M<sub>☉</sub> &gt; 25</td><td colspan="3">t/t<sub>MS</sub> ≤ 0.5</td><td colspan="3">t/t<sub>MS</sub> &gt; 0.5</td><td colspan="3"></td></tr><tr><td></td><td colspan="3">33.7</td><td colspan="3">33.6</td><td colspan="3">13</td></tr><tr><td></td><td colspan="3">⟨V sin i / V_c⟩</td><td colspan="3">⟨V sin i / V_c⟩</td><td colspan="3">⟨V sin i / V_c⟩</td></tr><tr><td>M/M<sub>☉</sub> ≤ 7.5</td><td>0.42</td><td>0.54</td><td></td><td>0.54</td><td>0.69</td><td></td><td>0.61</td><td>0.77</td><td></td></tr><tr><td>7.5 &lt; M/M<sub>☉</sub> ≤ 15</td><td>0.49</td><td>0.63</td><td></td><td>0.56</td><td>0.72</td><td></td><td>0.57</td><td>0.73</td><td></td></tr><tr><td>15 &lt; M/M<sub>☉</sub> ≤ 25</td><td>0.40</td><td>0.51</td><td></td><td>0.59</td><td>0.75</td><td></td><td>0.63</td><td>0.81</td><td></td></tr><tr><td>M/M<sub>☉</sub> &gt; 25</td><td colspan="3">t/t<sub>MS</sub> ≤ 0.5</td><td colspan="3">t/t<sub>MS</sub> &gt; 0.5</td><td colspan="3"></td></tr><tr><td></td><td colspan="3">0.56</td><td colspan="3">0.57</td><td colspan="3">0.72</td></tr></tbody></table>

**Notes.** For both series of average velocity ratios the standard deviations are ≥0.17.

inclinations, while $k < 0$ produces an excess of low inclination angles as compared to those predicted by randomly distributed inclined axes ($k = 0$). Figure 10a shows some functions $P_k(i)$, while Fig. 10b shows the histograms calculated for the same number of objects as in our stellar sample. The observed distribution $\Psi(v)$ of apparent ratios $v = V \sin i / V_c$ and the distribution $\Phi(u)$ of true ratios $u = V / V_c$ is now written as

$$ \Psi(v, k) = A_k v^{1+k} \int_v^{\infty} \Phi(u, k) u^{-(1+k)} (u^2 - v^2)^{-1/2} du. \quad (16) $$

Using the Lucy-Richardson algorithm, we obtain the distributions shown in Fig. 10c. In these calculations we used Stoeckley's corrections derived in Sect. 3.5.1. Although it is not shown explicitly, we note that the statistical uncertainty affecting the correction of distributions for measurement uncertainties is smaller than the spread of distributions $\Phi(u, k)$ obtained with $-0.5 \le k \le 0.5$. The modes and dispersions of functions $\Phi(u, k)$ due to $P_k(i)$ given in Eq. (15) are almost the same. On the contrary, their skewness and kurtosis are sensitive to $k$.

One might be tempted to inquire how the distribution of the inclination angles determined in this work compares with the random distribution. Using the derived inclinations and their uncertainties, we obtained the histogram shown in Fig. 11 (blue), which is compared with the histogram obtained for randomly distributed inclinations (white). In both cases the respective sampling uncertainties are shown. The student test performed on the individual classes in the interval $0 < \sin i \le 0.5$ shows that both distributions are the same to better than 99% confidence level, while they are definitely different in $0.5 \le \sin i \le 1.0$. In this last interval of inclinations, a function given by Eq. (15) with $k < 0$ might fit the distribution.

Before concluding anything about the actual distribution of the inferred angles $i$, it is important to note the following: (1) the $C_T$, $C_L$, and $C_G$ functions entering Eq. (3) do not vary strongly in the interval $0.7 \le \eta \le 1.0$, which according to the algorithm used to solve this system of equations, the propagation of the uncertainties affecting the observed parameters produces slightly overestimated values of $\eta$ and, thus, the number of stars with inclinations in the interval $35^\circ \le i \le 70^\circ$ is increased; (2) for some unknown reason, the stellar data used here do not frequently produce $\sin i \ge 0.95$, a high value indeed, but low enough to underestimate the number of stars with $i \ge 75^\circ$; (3) in spite of the care we payed to estimate the bolometric luminosity of stars, $L^{\text{app}}$, this parameter is probably the main source of uncertainties because we can never be absolutely sure that the effects due to the

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** a) Density probability distributions $P_k(i)$ against $\sin i$ for $k = -0.5, 0.0$ and $+0.5$. b) Frequency histograms of inclinations $i$ against $\sin i$ for $k = -0.5, 0.0$ and $+0.5$. c) Distributions $\Phi(u, k)$ of ratios $u = V / V_c$ assuming that the distributions of inclinations angles is given by $P_k(i)$ in Eq. (15). The tested cases are for $k = 0$ (random) and nonrandom distributions with $k = +0.5$ and $k = -0.5$.</figcaption></figure>

variable circumstellar envelope of Be stars have been correctly or entirely accounted for; and (4) all stars were assumed to be rigid rotators. Nevertheless, they can have some degree of differential rotation, which would then require that the line broadening be reinterpreted and given the right meaning to the $V \sin i$ parameter (see Paper II).

Thus, because of the difficulties listed above in determining the absolute inclination angles, it is not possible to figure out the actual distribution of inclinations of the studied Be stars. On the other hand, in Fig. 7 we have shown that the distribution of the velocity ratios $(V/V_c)_{pncr} = [(V \sin i)_{pncr} / V_c] / \sin i_{pncr}$ (crosses), which depends on the derived inclination angles $i$, fits very closely the distribution of ratios $V/V_c$ obtained assuming that the inclination angles are distributed at random. We can then conclude that even if the individual inclinations are still rough estimates, they are precise enough to warrant that the conclusions drawn from the statistical distribution of $V/V_c$ obtained in Sect. 3.5 is reliable.

# 4. Effect carried by the macroturbulence

## 4.1. Preliminary comments

Since the early fifties it has been recognized that spectral lines may undergo broadening due to random motions of eddies in
<figure><img src="image_12.png" /><figcaption>Fig. 11. Comparison of the distribution of inferred inclination angles for the program stars (cyan histogram) with the histogram for randomly distributed inclinations (white histogram). The error bars correspond to sampling uncertainties.</figcaption></figure>

the stellar atmospheres (Unsöld & Struve 1949; Huang & Struve 1953). When these movements are characterized by distance scales on the order or lower than the mean free-path of photons, they are called microturbulence. These movements affect the line absorption coefficient and produce the so-called “second-class” line broadening. Movements of eddies implying distance scales that are larger that photon mean free-paths are called macroturbulence and they do not change the effective mean atomic-line absorption coefficient. The line broadening is then produced by the Doppler effect associated with the macroscopic motion of eddies and is called first class. Howarth et al. (1997) claim that macroturbulence can be an important line broadening mechanism in O- and early B-type supergiants that adds to rotation. Nevertheless, its determination is not unique, since the same broadening can be produced by combining large ranges of macroturbulent velocities $v_{mt}$ and $V\sin i$ values (Ryan et al. 2002).

Like microturbulence, macroturbulence hides in principle a series of motions due to physical nature that is not yet well identified. Cantiello et al. (2009) evoke the subphotospheric convection to account for the microturbulent motions, while Aerts et al. (2009, 2014) suggest that macro-turbulence could be ascribed to low-amplitude gravity modes of nonradial pulsations. Similar conclusions are also put forward by Simón-Díaz (2015).

Two models are currently used to describe the line broadening carried by macroturbulent motions: isotropic Gaussian and anisotropic with radial and tangential components, which each have a Gaussian-dependent velocity distribution (Gray 1975, 1992). Although the anisotropic model depends on at least three free parameters, these authors assume that both models are characterized by a unique dispersion of macroturbulent velocities $\sigma_{mt}$, which is somewhat correlated with the stellar effective temperature. Generally, more effective fits of spectral line profiles are obtained using the anisotropic model, which requires slightly larger values of $v_{mt}$ than the isotropic model (Simón-Díaz & Herrero 2007, 2014; Dufton et al. 2006a). Both models are, however, nothing but spectral line fitting procedures that mask unidentified line broadening phenomena.

Up to now, macroturbulence was mostly explored in early-type stars of several luminosity classes, mainly supergiants, but all they have low values of $V\sin i$ (cf. Howarth et al. 1997; Ryan et al. 2002; Dufton et al. 2006a; Simón-Díaz & Herrero 2007; Fraser et al. 2010; Sundqvist et al. 2013). The existence of macroturbulence was also noticed in later B-type supergiants (Dufton et al. 2006a). According to Simón-Díaz & Herrero (2014) and Sundqvist et al. (2013) macroturbulence should not be entirely negligible in early-type dwarfs to giants. We can then expect that macroturbulent-like movements also exist in the atmospheres from B-type dwarfs to giants of lower effective temperatures than those explored in the above cited works. This can be the case of all those stars that undergo nonradial pulsations and/or a wide funnel of instabilities induced by the rotation. We assume then that the photospheric spectral lines currently used to determine the $V\sin i$ parameters of Be stars are affected by macroturbulence. It is worth noting that the macroturbulent motions are detected on average only when apparent rotational velocities are low, i.e., $V\sin i \le 150$ km s$^{-1}$. Spectral lines broadened by $V\sin i \ge 150$ km s$^{-1}$ remain fairly insensitive to broadening induced by macroturbulent motions. Although in these cases we are not able to detect such motions, we should not conclude that they are not present, regardless of their physical origin. This is particularly valid in Be stars, where atmospheres can undergo significant upheavals maintained by nonradial pulsations, a large spectrum of instabilities induced by the rapid rotation, and possible disordered magnetic fields generated by underphotospheric convection (Clement 1979; Maeder et al. 2008; Cantiello et al. 2009; Cantiello & Braithwaite 2011).

## 4.2. Formulation

Following the composition rules of probabilities given in Appendix D, the relation between the frequency $\Psi(v)$ of velocity ratios $v = (V\sin i)_{mt}/V_c$, where the $(V\sin i)_{mt}$ is the projected rotational velocity affected by macroturbulence, and the frequency $\Phi(u)$ of ratios $u = V\sin i/V_c$ corrected for macroturbulence, is given by

$$ \Psi(v) = \int_{0}^{v} \Phi[v - \epsilon(v_{mt}, u)] \left\{ \frac{\phi(v_{mt}, u)}{\partial\epsilon/\partial v_{mt}} \right\} du, \quad (17) $$

where we have written

$$ v = u + \epsilon(v_{mt}, u) \quad (18) $$

with $\epsilon(v_{mt}, u)$ representing the overestimation of the actual apparent rotation velocity assumed by the macroturbulence. In Eq. (17), $\phi(v_{mt})$ is the occurrence probability of the macroturbulent velocity $v_{mt}$. We assume that this distribution is the same for all B-type stars from dwarfs to giants.

The definition of $\epsilon(v_{mt}, u)$ given in Eq. (18) is the same as used by Simón-Díaz & Herrero (2007, see their Fig. 10), and it is adapted to account for the effect introduced by macroturbulence on the distribution function of rotational velocities. Its formal aspect can make us believe that a single scalar can also accurately account for the effects of macroturbulence on line profiles. In fact, the effects on a line profile depend on several parameters, depending on the approximation used to describe the macroturbulent motions: isotropic or anisotropic (Gray 1975, 1992). In any case, all these parameters are free without clear physical meaning and their values are not unique, as shown by Ryan et al. (2002). The parameter $\epsilon(v_{mt}, u)$ simply accounts for the difference between the $V\sin i$ parameters determined with the
<figure><img src="image_13.png" /><figcaption>Fig. 12. Normalized distributions $\Phi(v_{\text{mt}}, u)$ against the macroturbulent velocity $v_{\text{mt}}$ given by Eq. (19) for some $u = V\sin i/V_c$ values that are identified in the figure.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Macroturbulent broadening $\epsilon(v_{\text{mt}}, u)$ for 10 values of the macroturbulent velocity $v_{\text{mt}}$ scaled by steps of 10 km s$^{-1}$, as a function of the ratio $u = V\sin i/V_c$.</figcaption></figure>

FT method, once from rotationally broadened spectral lines convolved and then from not convolved line profiles for macroturbulence effects. We used the anisotropic approach to represent macroturbulence in the currently used monoparametric version (see Simón-Díaz & Herrero 2007). Other details for its use in Eq. (18) are given in Sect. 4.3.

To carry out the deconvolution of $\Phi(u)$ in Eq. (17) we need to specify $\phi(v_{\text{mt}})$. To this end we used the compilation of $v_{\text{mt}}$ values presented by Simón-Díaz (2015), which encompass 100 O- and 200 B-type stars of different spectral types and luminosity classes. This compilation is presented in Fig. 1 of Simón-Díaz (2015) as a two-dimensional distribution $\phi(v_{\text{mt}}, u)$ of discrete macroturbulent velocity $v_{\text{mt}}$ determinations as a function of $V\sin i$. Making cuts in this diagram by intervals of $V\sin i$ against $v_{\text{mt}}$, and cuts by intervals of $v_{\text{mt}}$ against the $V\sin i$, we produced least-squares fits of the resulting one-dimensional distributions of points with the following system of functions;

$$ 
\begin{aligned}
\phi(v_{\text{mt}}, u) &= [1 - \gamma(u)]\phi_G(v_{\text{mt}}) + \gamma(u)\phi_F(v_{\text{mt}}) \\
\phi_G(v_{\text{mt}}, u) &= c_g \exp\left\{-\frac{1}{2\sigma^2}[v_{\text{mt}} - v_M(u)]^2\right\} \\
\phi_F(v_{\text{mt}}, u) &= c_f \left\{1 - \exp\left[\frac{A(u)}{v_{\text{mt}}}\right]\right\}^{b(u)} \\
A(u) &= \left[\frac{b(u)}{b(u) - 1}\right]v_M(u),
\end{aligned}
\quad (19) $$

which we adopted as the analytical translation of the $\phi(v_{\text{mt}}, u)$ distribution of points in Simón-Díaz (2015). This density of points increases as a function of $v_{\text{mt}}$ and attains a maximum at $v_{\text{mt}} = v_M$, which in turn is an increasing function of $V\sin i$. For $V\sin i \gtrsim 80$ km s$^{-1}$ the points acquire a rather uniform distribution in the $(V\sin i, v_{\text{mt}})$ diagram. In Eq. (19), $c_g$ and $c_f$ are normalization constants. $\gamma(u)$ is a linear function of $u = V\sin i/V_c$ that varies from 0 to 1 as $V\sin i$ goes from 0 to 80 km s$^{-1}$. The exponent $b(u)$ accounts for the density of points, whose value ranges from $b=12$ in the interval $0 \le u \le 0.2$ ($0 \le V\sin i \le 80$ km s$^{-1}$) to about $b=100$ as $V\sin i \to 200$ km s$^{-1}$. The dispersion $\sigma = 17.5$ km s$^{-1}$ as well as the functions $v_M(u)$, $\gamma(u)$, and $b(u)$ are determined in the fitting procedure of the mentioned diagram in Simón-Díaz (2015). The obtained analytical distribution $\phi(v_{\text{mt}}, u)$ against $v_{\text{mt}}$ for some values $0.0 \le u = V\sin i/V_c \le 0.4$ is shown in Fig. 12.

### 4.3. Distribution of true velocities corrected from macroturbulent broadening

We calculated a great number of He I 4471 and Mg II 4481 line profiles for many effective temperatures and $\log g = 3.8$ to sketch an average MS evolutionary phase. They were all broadened assuming that stars were gravity-darkened rigid rotators (Frémat et al. 2005). The obtained line profiles were then broadened for macroturbulence with the currently used simplified anisotropic model. For each adopted macroturbulent velocity $v_{\text{mt}}$, we calculated the weighted average of the broadening differences $\epsilon(v_{\text{mt}}, u)$. The weights are the fractions of sample stars that enter per interval of effective temperatures. Figure 13 shows the deviations $\epsilon(v_{\text{mt}}, u)$ obtained for several values of $v_{\text{mt}}$ against $u = V\sin i/V_c$.

To obtain $\Phi(u)$ in Eq. (17), we used the Lucy-Richardson algorithm (Lucy 1974; Richardson 1972) and the probability distribution $\phi(v_{\text{mt}}, u)$ given by Eq. (19). The results thus obtained are shown in Fig. 14. The blue line represents the normalized distribution of parameters $V\sin i/V_c$ corrected for both, measurement uncertainties and Stoeckley's effect (cf. Fig. 2). The red curve represents the normalized distribution $\Phi(u)$ of $u = V\sin i/V_c$ statistically corrected for macroturbulent line broadening obtained with $b=12$ in $0 \le u \le 0.2$. As expected, the distribution is slightly affected in the interval $0.0 \le u = V\sin i/V_c \le 0.4$, and suggests a modest increase of genuine low rotators, while the mode of the distribution is slightly reduced. We do not have a reliable explanation for the bump-shaped excess at $u = V\sin i/V_c \sim 0.10-0.15$. If there were a higher number of points in $0 \le u \le 0.2$ than actually displayed by Simón-Díaz (2015), we should need $b < 12$. Figure 14 shows the effect produced by a hypothetical exponent $b=6$ (black dashed curve).

In conclusion, owing to the correction for macroturbulent line broadening, the number of low rotators is slightly increased to the detriment of rotators with velocity ratios near the mode of the distribution, but the mode itself is almost preserved.

## 5. Effect carried by the presence of a binary component

In Be stars, binarity not only helps to form a circumstellar disk during a mass-transfer episode (Kriz & Harmanec 1975; Harmanec 1987), but also to spin up one of the components owing to the angular momentum gained during a mass transfer
<figure><img src="image_15.png" /><figcaption>Fig. 14. Distribution $\Phi(u)$ after correction for macroturbulence. The blue line represents the normalized distribution $\Psi(v)$ of parameters $v = V\sin i/V_c$ corrected for both, measurement uncertainties and Stoeckley’s effect. The red and dashed black curves represent the normalized distribution $\Phi(u)$ of ratios $u = V\sin i/V_c$ corrected for macroturbulent motions using different values of $b$ (see text).</figcaption></figure>

event (Packet 1981; Gies 2000; Harmanec et al. 2002; de Mink et al. 2013). In the context of our subject, it has long been accepted that binarity can alter the $V\sin i$ parameters of the individual components.

As in other stellar classes, a rather large fraction $f_{\text{bin}}$ of Be stars studied in this work can be binaries. Recent statistics show that among massive stars this frequency can be as high as 70% (de Mink et al. 2013), while it climbs to nearly 100% (Sana et al. 2012, 2014) for O-type stars. Chini et al. (2012) report that the frequency of binaries is higher than 80% for masses $M > 16 M_{\odot}$, and it steadily decreases to 20% for masses $M \sim 3 M_{\odot}$. Oudmaijer & Parr (2010) find that the frequency of binaries among B stars without emission is 30% and it is the same for Be stars. In the past there have been several surveys of B star binarity: Abt et al. (1990) obtain 74% binaries among B2-B5-type stars, considering both spectroscopic and visual; Kouwenhoven et al. (2007) infer that in the Sco OB2 association there must be 70% binaries among intermediate-mass stars (B- to A-type); and Mason et al. (2009) find 64% binaries in a speckle interferometry survey of Galactic B stars.

Beyond the low-mass extreme of B-type stars, Fuhrmann & Chini (2012) report that the percentage of binary stars may probably increase to 100% from F to late A-type stars, although these last are considered with all their aspects of multiplicity. Recent studies based on adaptive optics surveys (De Rosa et al. 2014), and VLTI/PIONIER interferometry (Marion et al. 2014) note an increase of 13% A-type stars among objects considered single object up to now.

Since the frequency of binary systems for O-type stars approaches 100% and the multiplicity of all kinds among late A-type stars attains similar percentages, we can ask whether the above noted low frequency of binary systems of B spectral type is biased in some way.

According to environments, authors and methods used to detect binaries, the frequency for B-type binaries ranges from 30% to 70%. Frequently, the search of binaries among them is

limited to the use of Balmer lines, whose large broadening easily hides the existence of components leading thus to underestimate the frequency of binaries. On the other hand, the occurrence of the Be phenomenon is the largest among stars of apparent B2 spectral type. These stars have masses of apparent O9-B1 spectral types (Zorec & Briot 1997; Zorec et al. 2007), implying $M \ge 12 M_{\odot}$, which enters the domain of massive stars, where binarity is frequent. As a consequence of these uncertainties, the curious lack of binaries among B-type stars appears rather doubtful. Thus, we proceed by making the ansatz that at least 70% of Be stars are binaries.

## 5.1. Formulation

We assume that the measured $V\sin i$ parameter of an object in a binary system mainly reflects the rotation of the primary. This value can then be affected by the orbital angular velocity, but the effect depends on the mass ratio and the separation of components.

The radial velocity of a point $\bar{P}$ located on the observed hemisphere of the tested primary component, whose rotation axis is for simplicity taken parallel to that of the orbital rotation, is given by (Limber 1963; Moreno & Koenigsberger 1999; Moreno et al. 2005)

$$V_x = \left\{ \frac{q}{1+q} r(\theta)\omega(t) \left[ r(\theta) \frac{e}{p} r(\theta) \sin^2\theta + \cos\theta \right] - \left[ \Omega_p + \omega(t) \right] x_1 \right\} \sin i, \quad (20)$$

where $q = M_s/M_p$ is the mass ratio of the secondary to the primary; $\omega(t)$ is the time dependent orbital angular velocity; $e$ is the eccentricity of the orbit; $p = a(1 - e^2)$ is the conic parameter with $a = a_p + a_s$ the semimajor axis of the relative orbit, with $a_p$ and $a_s$ being those of the components; $\theta$ is the polar azimuthal coordinate; $\Omega_p$ is the value of the nonsynchronous angular velocity of the primary; $r(\theta)$ is the time-dependent separation of the components; $x_1$ is the coordinate of the test point $\bar{P}$ measured positively from the center of the primary toward the secondary on the axis that joins both components; and $\sin i$ accounts for the inclination between the rotation axis and the direction toward the observer. Because of several other uncertainties that we may deal with in the present discussion, we assume that all binaries have circular orbits ($e = 0$), so that $\omega(t) = \omega = \text{constant}$. We also do not take into account the periodical apparent changes of $V\sin i$ due to spectral effects produced by mutual reflection effects. The determination of $V\sin i$ requires that all spectral lines are reduced to a common radial velocity. This means that the first term in Eq. (20), which represents the orbital motion, is removed. According to these simplifications, the effective velocity that broadens the spectra lines reduces to (cf. Eq. (20))

$$V_x = [\Omega_p + \omega] x_1 \sin i, \quad (21)$$

which means that $V$ in the measured $V\sin i$ is

$$V = (\Omega_p + \omega)R_p = V_p + \Delta V, \quad (22)$$

where $R_p$ is the radius of the primary. Using the third law of Kepler, we readily obtain the expression for the normalized perturbation component $\Delta v = \Delta V/V_c$,

$$\Delta v = [(1 + q)\delta^{-3}]^{1/2}, \quad (23)$$
<figure><img src="image_16.png" /><figcaption>Fig. 15. Curves log f(Δv). All 9 solutions obtained with -0.6 ≤ α ≤ 0.3 and -0.5 ≤ β ≤ 0.5 in Eq. (26) for f(Δv) are comprised in the gray-shaded region. The red line represents the average of the 9 tested solutions and the blue line corresponds to the solution for α = β = 0.</figcaption></figure>

where we have written $\delta = a/R_p$. It is easy to see that the range of variation of Δv is limited to 0 ≤ Δv ≤ 0.5, where Δv = 0.5 corresponds to the extreme configuration given by q = 1 and $\delta_{min} = 2$.

The calculation of the effect carried by Δv on the observed distribution of V/Vc depends on the occurrence probability f(Δv) of Δv. Following the definition of probability and the relations presented in Appendix D, the frequency function f(Δv) can be obtained with the surface integral

$$f(\Delta v)d(\Delta v) = \iint_{\Delta\sigma} Q(q)D(\delta) dq d\delta, \qquad (24)$$

which is calculated over the area Δσ and limited by the curves Δv = Δv(q, δ) and Δv + d(Δv) = Δv(q, δ). In Eq. (24), Q(q) is the probability distribution of the binary mass-ratios q, while D(δ) represents the probability of the occurrence of the orbital separation δ. We adopt for Q(q) and D(δ) the following expressions frequently used in the literature (Kobulnicky & Fryer 2007; Reggiani & Meyer 2013; de Mink et al. 2013),

$$\begin{cases} Q(q) = q^{\alpha} \\ D(\delta) = \frac{(\ln \delta)^{\beta}}{\delta} \end{cases} \qquad (25)$$

Introducing these definitions into Eq. (24) and knowing that for each δ = const. it is dq = d(Δv)/|Δv/Δq|, where 1/(Δv/Δq) = 2Δvδ³, the integration of Eq. (24) over δ from Δv⁻²/³ to (2Δv²)¹/² (see Eq. (23)), we obtain

$$f(\Delta v) = \frac{f_0}{\Delta v} \int_0^1 q^{\alpha} \left[ \ln \left( \frac{1+q}{\Delta v^2} \right) \right]^{\beta} dq, \qquad (26)$$

where f₀ stands for a normalization constant. Since the exponents α and β are highly uncertain (Shatsky & Tokovinin 2002; Kobulnicky & Fryer 2007; Reggiani & Meyer 2013; de Mink et al. 2013, and references therein) to account for the large range of their possible values, we calculated Eq. (26) using nine pairs (α, β) with α = -0.6, -0.15, 0.3 and β = -0.5, 0.0, 0.5. The nine functions f(Δv) thus obtained are all comprised in the shaded region shown in Fig.15. Because the higher frequencies f(Δv) are for 0 ≤ Δv ≤ 0.2, we can safely adopt a unique average relation,

<figure><img src="image_17.png" /><figcaption>Fig. 16. Distribution functions of true velocity ratios V/Vc. Ψ(v)Σ+mt is the distribution of v = V/Vc corrected for measurement uncertainties, Stoeckley's effect (Σ), macroturbulent velocities (mt) with b = 12 and projection factor sin i. It distinguishes thus from Φ(u)b=12 shown in Fig. 14, which is not corrected for effects owing to sin i. Full red line corresponds to the distribution of u = V/Vc corrected for Σ + mt and binarity with f_bin = 0.7. Dashed red line corresponds to the distribution corrected for Σ + mt and binarity assuming f_bin = 1.0.</figcaption></figure>

which nearly coincides with the expected canonical solution for α = β = 0 (blue curve in Fig.15) adopted by de Mink et al. (2013), which can be nicely represented by the following simple interpolation expression:

$$f(\Delta v) \approx \frac{0.168}{\Delta v}. \qquad (27)$$

Since Q(q) and D(δ) are given by power laws, the normalization of f(Δv) must be carried out with limited values of Δv. Apart from Δvmax = 0.5, the other limiting value Δv → 0 was here assumed to be provided by Δvmin = 10⁻⁴, which implies a negligible effect on V and roughly corresponds to an orbital separation δmax = 600.

## 5.2. Distribution velocities corrected for effects due to binarity

Assuming that the same fraction f_bin of binaries occurs in all spectral types and luminosity classes, the relation between the observed frequency function Ψ(v) of ratios v = V/Vc perturbed by the presence of binaries and the “actual” frequency Φ(u) of “unperturbed” ratios u = v - Δv is given by

$$\Psi(v) = f_{\text{bin}} \int_{0}^{v} \Phi(u) f(v-u) du + (1 - f_{\text{bin}}) \Phi(v), \qquad (28)$$

which is solved with the Lucy-Richardson deconvolution algorithm (Richardson 1972; Lucy 1974). Figure 16 shows the results obtained for f_bin = 0.7 and f_bin = 1.0. In this figure, Ψ(v)Σ+mt stands for the input distribution of velocity ratios V/Vc corrected for Stoeckley's effect and macroturbulent velocities. Apart from the modes, which are systematically shifted to lower values the higher is the fraction f_bin, there is a small secondary bump at V/Vc ~ 0.1 that appears due to the correction for macroturbulent velocities. Perhaps, according to speculations in Fig. 14,
instead of a “bump” at $V/V_c \sim 0.1$, there is a “depression” in $0.15 \le V/V_c \le 0.35$, which reveals a possible insufficient correction for macroturbulence or even other perturbations that we are not able to detect at such rotational velocity ratios.

## 6. Comments and conclusions

The main purpose of this work is to obtain the distribution $\Phi(u)$ of rotational velocity ratios $u = V/V_c$ ($V_c$ is the critical velocity at the equator) of a sample of Galactic classical Be stars. The distributions of ratios $V/V_c$ were corrected statistically for a number of phenomena affecting the value of the $V\sin i$ that blurs the information on the stellar rotation contained in the factor $V$. The function $\Phi(u)$ obtained in the present work is used in the second paper of this series (Zorec et al. 2016) to inquire about the average properties of the surface rotation in Be stars.

Because the parameter $V\sin i$ of Be stars suffers from effects of gravitational darkening (GD), we discussed two statistical deconvolution procedures that take “Stoeckley’s underestimation” in different ways and produce distributions of ratios $V/V_c$ from the observed distributions of projected velocity ratios $V\sin i/V_c$. The results obtained using these methods can be summarized as follows:

1. We adopted a sample of 233 classic Be stars that was previously studied in Chauville et al. (2001), Frémat et al. (2005), Zorec et al. (2005), Frémat et al. (2006) and Levenhagen & Leister (2006), whose apparent (viewing angle dependent) fundamental parameters were rediscovered. The selection of this stellar sample is justified by the similarity in the approaches and in precautions taken to determine the $V\sin i$ (Sect. 2.1).

2. The translation of the distribution of $V\sin i/V_c$ into another one corrected for the $i$-viewing angle effect proceeds by assuming that $i$ is distributed at random and considering that the parameters $V\sin i$ are systematically underestimated because of GD (Sect. 3).

As it is still rather widely presumed in the literature that Be stars are critical rotators, we used the original von Zeipel formulation to determine the distribution $\Phi(V/V_c)$ characterized by the largest possible mode to calculate the GD effect, i.e., the larger value of the ratio $V/V_c$ at which $\Phi$ takes its maximum.

3. We discussed two deconvolution methods:

(a) The Cranmer (2005) method (see Sect. 3.4) uses a Monte Carlo simulation based on a mathematical stellar sample, where each parameter $V\sin i$ is corrected for the GD effect, and sampled statistically according to the quantities $V$ and $i$. This method also uses an imposed analytical expression for $\bar{\Phi}$ to represent the actual distribution $\Phi$ whose defining free parameters are selected by searching the best fit of the observed distribution $\Psi$ of $V\sin i$. We have shown that even when the set of moments that characterizes each function $\bar{\Phi}$ is similar to that of the actual distribution $\Phi$, the Kolmogorov-Smirnov test indicates that not all analytical functions $\bar{\Phi}$ defined with the same number of free parameters can be used to represent the sought function $\Phi$ reliably. In this context, the trapezoidal function $\bar{\Phi}$ is less appropriate. From this study it also emerges that there can be analytical expressions for $\bar{\Phi}$ where the characterizing free parameters can have different values, but they lead to equally “well satisfying fits” of the observed distribution $\Psi$.

(b) The Lucy-Richardson (1972, 1974) method (see Sect. 3.5) is based on the iteration of the distribution $\Phi(u)$ of ratios $u = V/V_c$ and the distribution $\Psi(v)$ of ratios $v = V\sin i/V_c$. This method is free from any imposed analytical expression of $\Phi$, but requires that the ratios $V\sin i/V_c$ be corrected for Stoeckley’s underestimation in advance. Such a correction requires a thorough determination of the fundamental parameters of each studied object. To this end, we introduced a method to determine the stellar mass, age, rotational ratio $\eta$ (ratio of centrifugal to the gravitational acceleration at the equator), and approximate inclination angle $i$ of each star in the studied sample (Sect. 3.5.1). The distribution $\Phi(u)$ of ratios $u = V/V_c$ corrected for GD effect finally adopted here is that obtained with the Lucy-Richardson method.

4. We separated the stars into three MS evolutionary groups, early, intermediate, and late, to provide an idea on the evolution effect on the distribution of ratios $V/V_c$ (Sect. 3.6). In spite of the uncertainties related to the definition of reliable statistical samples, we found that once the mode ($V/V_c \lesssim 0.77$) is attained by Be stars at some intermediate evolutionary phase in the MS, it will stay unchanged up to the end of the MS evolutionary lifespan. This behavior may represent a property that can be specific to Be stars.

5. We attempted a discussion of the effect introduced on $\Phi(V/V_c)$ by a nonrandom distribution of the inclination angle $i$ of the stellar rotation axis (Sect. 3.7). According to the assumed function that gives the probability of occurrence of inclination angles, the mode of $\Phi$ does not change heavily but its skewness and kurtosis depend on the parameter that represents the deviation to the randomness.

The Lucy-Richardson iteration method was then used to introduce two different corrections to the distributions of velocity ratios to account for phenomena that carry systematic overestimations of parameters $V\sin i$:

I. We have taken into account the effects on the distribution $\Psi(V\sin i/V_c)$ due to the macroturbulent motions in the stellar atmospheres (Sect. 4). This correction slightly increases the frequency $\Phi(V/V_c)$ in the interval of velocity ratios $0 \le V/V_c \le 0.3$, but does not significantly change the mode of $\Phi(V/V_c)$.

II. As there can be an additional component in the $V\sin i$ owing to the orbital motion of objects in a binary system (Sect. 5), we formulated the correction for this effect and considered that there must be at least $f_{\text{bin}} = 70\%$ binary stars among Be stars. This $f_{\text{bin}}$ is because there are no clear reasons that among B-type stars the frequency of binaries is significantly smaller than in both B-type spectral extremes, i.e., more and less massive objects than those assigned to B-type objects. The effect for binarity is the last correction introduced here on the distribution of velocity ratios $V/V_c$.

Thus, according to the corrections taken into account thus far, we conclude the following:

(i) The mode of $\Phi(V/V_c)$ depends on the frequency $f_{\text{bin}}$ of binaries: for $f_{\text{bin}} = 0.7$ the mode is at $(V/V_c)_{\text{mode}} = 0.663$, while for $f_{\text{bin}} = 1.0$ it becomes $(V/V_c)_{\text{mode}} = 0.644$. If we had admitted that $f_{\text{bin}} = 0.3$, the mode would be at $(V/V_c)_{\text{mode}} = 0.680$, while for $f_{\text{bin}} = 0.0$ it is at $(V/V_c)_{\text{mode}} = 0.681$. At fractions $f_{\text{bin}} \ge 0.7$ the number of rotators with $0.75 \le V/V_c \le 1.0$ reduces significantly.
(ii) In Fig. 16 we can see that for $f_{\text{bin}} = 0.7$ the most frequent rotation of classic Be stars is characterized by $(V/V_c)_{\text{mode}} = 0.663$, which implies $\eta = 0.35$ or $\Omega/\Omega_c = 0.84$. From $\Phi(u)_{\text{\Sigma+mt+bin0.7}}$ we can show that 90% of program stars, i.e., 210 objects, rotate with $0.44 \le V/V_c \le 0.89$, which means $0.61 \le \Omega/\Omega_c = 0.98$ or $0.14 \le \eta \le 0.73$. Accordingly, 7% of them (16 stars) have $V/V_c \le 0.4$ and barely 3% of them (7 stars) have $V/V_c \ge 0.9$. Since half of the Be stars in the whole sample rotate roughly with $V/V_c \le 0.65$, saying that Be stars are critical rotators sounds excessive.

(iii) According to Fig. 9, the larger fraction of Be stars with $0.44 \le V/V_c \le 0.89$ is in the second half of the MS evolutionary phase ($t/t_{\text{MS}} \ge 0.5$).

(iv) From Fig. 16 we can see that the difference in velocities owing to the presence of binary components represent $\Delta V \approx 14 \text{ km s}^{-1}$ at $V/V_c = 0.4$ and $\Delta V \approx 38 \text{ km s}^{-1}$ at $V/V_c = 0.9$, which are measurable.

We note that Eq. (8) is not appropriate even for rigid rotators because $\beta_1$ is a function of the stellar colatitude $\theta$ and its value decreases as $\eta \to 1$ (Espinosa Lara & Rieutord 2011; Zorec et al. 2016). However, in the present contribution Stoeckley's correction $\Sigma(\eta, i, M, t)$ given in Eq. (7) was calculated imposing $\beta_1 = 1$, so as to overestimate the GD effect on purpose. This means that the correction for GD used in this work is the highest we could apply to an observed $V\sin i$ parameter and that the characteristic $V/V_c$ ratios of Be stars must actually be lower than inferred in this work.

Finally, if critical rotation is a transient phenomenon that helps the star to eject the matter needed to form its CE, the transition from some $V/V_c = 0.9$ to $V/V_c = 1.0$ would require that during a short time a sufficient amount of energy be provided to photospheric layers to let them at critical rotation. Critical rotation maintains, however, the surface layers bound to the star, which can be freed only if the rotation rate attains the escape velocity, $V_{\text{esc}} = \sqrt{2V_c}$.

Assuming that the excess of energy necessary to attain the critical rotation is supplied by a magnetic field, in a star with an equatorial rotation velocity $V$, the strength of this field must obey $B^2/2\mu_o \sim \rho_{\text{phot}}[1-(V/V_c)^2]/2$, where $\rho_{\text{phot}}$ is an average density of photospheric layers of dwarf B-type stars. In a typical Be star $V_c = 480 \text{ km s}^{-1}$ and $\rho_{\text{phot}} \sim 6 \times 10^{-10} \text{ g cm}^{-3}$. The above mentioned change of surface rotational velocities would then imply $B \ge 2 \times 10^3 \text{ G}$. Such magnetic fields could certainly be detected if they existed and were global. As they are not detected, magnetic fields of such intensities could nonetheless exist, but must be concentrated in small areas, and maintained perhaps by a dynamo process hosted in the rotationally enlarged subsurface convection zones, such as those suggested by Clement (1979) and Maeder et al. (2008). Under such cases it would also be possible to think of activities that evoke flares with huge mass-ejection episodes through magnetic recombination phenomena. For example, emerging magnetic fields from the surface of some $10^2 \text{ G}$ by magnetic buoyancy in massive stars that originate in the thin “iron convection zone” were already foreseen by Cantiello & Braithwaite (2011).

In Paper II we complete the discussion started in the present work. To this end we assume that stars undergo surface differential rotation and derive: (1) the distribution function $\Phi(u)$ corrected for differential rotation; (2) the density probability function for the occurrence of a gradient of the angular velocity in the stellar surface that is consistent with the predicted line broadening and with the distribution $\Phi(u)$ of ratios $u = V/V_c$ obtained in Sect. 5.1. In Paper II we also discuss a number of measurement

and conceptual uncertainties affecting the $V\sin i$ determination of Be stars. These include the effects due to surface differential rotation, bivalued relations between line broadening and the $V\sin i$, expansion velocities in atmospheres of active stars, and effects on spectral lines due to tides in binaries.

*Acknowledgements.* This research has widely made use of the SIMBAD database, operated at CDS, Strasbourg, France. We are grateful to an anonymous referee for his/her detailed reading of our manuscript and for his/her helpful comments and suggestions.

## References

Abt, H. A., Gomez, A. E., & Levy, S. G. 1990, ApJS, 74, 551

Adelman, S. J., Pintado, O. I., Nieva, M. F., Rayle, K. E., & Sanders, Jr., S. E. 2002, A&A, 392, 1031

Aerts, C., Puls, J., Godart, M., & Dupret, M.-A. 2009, A&A, 508, 409

Aerts, C., Simón-Díaz, S., Groot, P. J., & Degroote, P. 2014, A&A, 569, A118

Aidelman, Y., Cidale, L. S., Zorec, J., & Arias, M. L. 2012, A&A, 544, A64

Aller, L. H., Appenzeller, I., Baschek, B., et al. 1982, Landolt-Börnstein: Numerical Data and Functional Relationships in Science and Technology – New Series, Gruppe/Group 6 Astronomy and Astrophysics, Vol. 2 Schaifers/Voigt: Astronomy and Astrophysics, Stars and Star Clusters

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Ballereau, D., Chauville, J., & Zorec, J. 1995, A&AS, 111, 423

Balona, L. A. 1975, MNRAS, 173, 449

Beeckmans, F., & Hubert-Delplace, A. M. 1980, A&A, 86, 72

Bodenheimer, P. 1971, ApJ, 167, 153

Bowman, A. W., & Azzalini, A. 1997, Applied smoothing techniques for data analysis: the kernel approach with S-plus illustrations, Oxford Stat. Sci. Ser. 18 (Clarendon Press)

Briot, D., & Zorec, J. 1987, in Physics of Be Stars, eds. A. Slettebak, & T. P. Snow, IAU Colloq., 92, 214

Brown, J. C., & Wood, K. 1992, A&A, 265, 663

Cameron, C., Saio, H., Kuschnig, R., et al. 2008, ApJ, 685, 489

Cantiello, M., & Braithwaite, J. 2011, A&A, 534, A140

Cantiello, M., Langer, N., Brott, I., et al. 2009, A&A, 499, 279

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1988, ApJ, 329, L33

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Chalonge, D., & Divan, L. 1952, Ann. Astrophys., 15, 201

Chandrasekhar, S., & Münch, G. 1950, ApJ, 111, 142

Chauville, J., Zorec, J., Ballereau, D., et al. 2001, A&A, 378, 861

Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925

Cidale, L., Zorec, J., & Tringaniello, L. 2001, A&A, 368, 160

Claret, A. 2012, A&A, 538, A3

Clement, M. J. 1979, ApJ, 230, 230

Cochetti, Y. R. 2014, Tésis de Licenciatura, Universidad Nacional de La Plata, Argentina

Collins, II, G. W. 1987, in Physics of Be Stars, eds. A. Slettebak, & T. P. Snow, IAU Colloq., 92, 3

Collins, II, G. W., & Truax, R. J. 1995, ApJ, 439, 860

Collins, II, G. W., Truax, R. J., & Cranmer, S. R. 1991, ApJS, 77, 541

Cranmer, S. R. 2005, ApJ, 634, 585

Dall, T. H., & Sbordone, L. 2011, J. Phys. Conf. Ser., 328, 012016

de Jager, C., & Nieuwenhuijzen, H. 1987, A&A, 177, 217

de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764, 166

de Rosa, R. J., Patience, J., Wilson, P. A., et al. 2014, MNRAS, 437, 1216

de Wit, W. J., Lamers, H. J. G. L. M., Marquette, J. B., & Beaulieu, J. P. 2006, A&A, 456, 1027

Deutschman, W. A., Davis, R. J., & Schild, R. E. 1976, ApJS, 30, 97

Divan, L., & Zorec, J. 1982a, in The Scientific Aspects of the Hipparcos Space Astrometry Mission, eds. M. A. C. Perryman, & T. D. Guyenne, ESA SP, 177, 101

Divan, L., & Zorec, J. 1982b, in Be Stars, eds. M. Jaschek, & H.-G. Groth, IAU Symp., 98, 61

Divan, L., Zorec, J., & Briot, D. 1982, in IAU Symp. 98, Be Stars, eds. M. Jaschek, & H.-G. Groth, 53

Divan, L., Zorec, J., & Andrillat, Y. 1983, A&A, 126, L8

Domiciano de Souza, A., Kervella, P., Moser Faes, D., et al. 2014, A&A, 569, A10

Dufton, P. L., Ryans, R. S. I., Simón-Díaz, S., Trundle, C., & Lennon, D. J. 2006a, A&A, 451, 603

Dufton, P. L., Smartt, S. J., Lee, J. K., et al. 2006b, A&A, 457, 265

Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146
Espinosa Lara, F., & Rieutord, M. 2011, A&A, 533, A43

Fitzpatrick, E. L., & Massa, D. 2005, AJ, 129, 1642

Floquet, M., Hubert, A. M., Hirata, R., et al. 2000, A&A, 362, 1020

Flower, P. J. 1996, ApJ, 469, 355

Fraser, M., Dufton, P. L., Hunter, I., & Ryans, R. S. I. 2010, MNRAS, 404, 1306

Frémat, Y., Zorec, J., Hubert, A.-M., & Floquet, M. 2005, A&A, 440, 305

Frémat, Y., Neiner, C., Hubert, A.-M., et al. 2006, A&A, 451, 1053

Friedemann, C., & Roeder, U.-K. 1987, Astron. Nachr., 308, 41

Friedemann, C., Guertler, J., Schielicke, R., & Dorschner, J. 1983, Astron. Nachr., 304, 237

Fuhrmann, K., & Chini, R. 2012, ApJS, 203, 30

Georgy, C., Granada, A., Ekström, S., et al. 2014, A&A, 566, A21

Gies, D. R. 2000, in The Be Phenomenon in Early-Type Stars, eds. M. A. Smith, H. F. Henrichs, & J. Fabregat, IAU Colloq., 175, ASP Conf. Ser., 214, 668

Gillich, A., Deupree, R. G., Lovekin, C. C., Short, C. I., & Toqué, N. 2008, ApJ, 683, 441

Gkouvelis, L., Fabregat, J., Zorec, J., et al. 2016, A&A, 591, A140

Gray, D. F. 1975, ApJ, 202, 148

Gray, D. F. 1992, The observation and analysis of stellar photospheres Camb. Astrophys. Ser. (Cambridge University Press)

Guertler, J., Schielicke, R., Dorschner, J., & Friedemann, C. 1982, Astron. Nachr., 303, 105

Hadrava, P. 1992, A&A, 256, 519

Hanuscihk, R. W., Dachs, J., Baudzus, M., & Thimm, G. 1993, A&A, 274, 356

Harmanec, P. 1987, in Physics of Be Stars, eds. A. Slettebak, & T. P. Snow, IAU Colloq. 92, 339

Harmanec, P., Bisikalo, D. V., Boyarchuk, A. A., & Kuznetsov, O. A. 2002, A&A, 396, 937

Howarth, I. D., Siebert, K. W., Hussain, G. A. J., & Prinja, R. K. 1997, MNRAS, 284, 265

Huang, S.-S., & Struve, O. 1953, ApJ, 118, 463

Huang, W., & Gies, D. R. 2006, ApJ, 648, 580

Hubert, A. M., Floquet, M., & Zorec, J. 2000, in The Be Phenomenon in Early-Type Stars, eds. M. A. Smith, H. F. Henrichs, & J. Fabregat, IAU Colloq., 175, ASP Conf. Ser., 214, 348

Hunter, I., Lennon, D. J., Dufton, P. L., et al. 2008, A&A, 479, 541

Jamar, C., Macau-Hercot, D., Monfils, A., et al. 1976, Ultraviolet bright-star spectrophotometric catalogue, A compilation of absolute spectrophotometric data obtained with the Sky Survey Telescope (S2/68) on the European Astronomical Satellite TD-1 (Paris: European Space Agency)

Jaschek, M., Slettebak, A., & Jaschek, C. 1981, Be Star Newsletter, 4, 9

Johnson, H. L., & Mitchell, R. I. 1975, Rev. Mex. Astron. Astrofis., 1, 299

Keller, S. C., Bessell, M. S., Cook, K. H., Geha, M., & Syphers, D. 2002, AJ, 124, 2039

Kippenhahn, R. 1977, A&A, 58, 267

Kobulnicky, H. A., & Fryer, C. L. 2007, ApJ, 670, 747

Kouwenhoven, M. B. N., Brown, A. G. A., Portegies Zwart, S. F., & Kaper, L. 2007, A&A, 474, 77

Kriz, S., & Harmanec, P. 1975, Bull. Astron. Inst. Czechoslovakia, 26, 65

Kuiper, G. P. 1935, PASP, 47, 15

Kurucz, R. 1993, Opacities for Stellar Atmospheres: [+0.0], [+0.5], [+1.0]. Kurucz CD-ROM No. 2 (Cambridge, Mass.: Smithsonian Astrophysical Observatory), 2

Lamers, H. J. G. L. M., Zickgraf, F.-J., de Winter, D., Houziaux, L., & Zorec, J. 1998, A&A, 340, 117

Lang, K. R. 1992, Astrophysical Data I. Planets and Stars (Springer-Verlag)

Levenhagen, R. S., & Leister, N. V. 2006, MNRAS, 371, 252

Limber, D. N. 1963, ApJ, 138, 1112

Lovekin, C. C., Deupree, R. G., & Short, C. I. 2006, ApJ, 643, 460

Lucy, L. B. 1967, Z. Astrophys., 65, 89

Lucy, L. B. 1974, AJ, 79, 745

Macau-Hercot, D., Jamar, C., Monfils, A., et al. 1978, Supplement to the ultraviolet bright-star spectrophotometric catalogue, A compilation of absolute spectrophotometric data obtained with the sky survey telescope (S2/68) on the European Astronomical Satellite TD-1 (Paris: European Space Agency), ESASR

Maeder, A. 1999, A&A, 347, 185

Maeder, A. 2009, Phys., Formation and Evolution of Rotating Stars (Springer)

Maeder, A., Georgy, C., & Meynet, G. 2008, A&A, 479, L37

Marion, L., Absil, O., Ertel, S., et al. 2014, A&A, 570, A127

Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009, AJ, 137, 3358

Meilland, A., Stee, P., Zorec, J., & Kanaan, S. 2006, A&A, 455, 953

Meilland, A., Millour, F., Kanaan, S., et al. 2012, A&A, 538, A110

Mennickent, R. E., Pietrzyński, G., Gieren, W., & Szewczyk, O. 2002, A&A, 393, 887

Milne, E. A. 1923, MNRAS, 83, 118

Mokiem, M. R., de Koter, A., Evans, C. J., et al. 2006, A&A, 456, 1131

Moon, T. T., & Dworetsky, M. M. 1985, MNRAS, 217, 305

Moreno, E., & Koenigsberger, G. 1999, Rev. Mex. Astron. Astrofis., 35, 157

Moreno, E., Koenigsberger, G., & Toledano, O. 2005, A&A, 437, 641

Moujtahid, A. 1998, Ph.D. Thesis, Université Pierre & Marie Curie, Paris

Moujtahid, A., Zorec, J., & Hubert, A. M. 1999, A&A, 349, 151

Moujtahid, A., Zorec, J., Hubert, A. M., Garcia, A., & Burki, G. 1998, A&AS, 129, 289

Moujtahid, A., Zorec, J., & Hubert, A. M. 2000a, in The Be Phenomenon in Early-Type Stars, eds. M. A. Smith, H. F. Henrichs, & J. Fabregat, IAU Colloq. 175, ASP Conf. Ser., 214, 506

Moujtahid, A., Zorec, J., & Hubert, A. M. 2000b, in I 175: The Be Phenomenon in Early-Type Stars, eds. M. A. Smith, H. F. Henrichs, & J. Fabregat, IAU Colloq. 175, ASP Conf. Ser., 214, 510

Napiwotzki, R., Schoenberner, D., & Wenske, V. 1993, A&A, 268, 653

Nasseri, A., Dembsky, T., Drass, H., Hoffmeister V. H., & Chini R. 2013, Central European Astrophysical Bulletin, 37, 51

Negruelua, I. 2004, Astron. Nachr., 325, 380

Nieva, M.-F. 2013, A&A, 550, A26

Oudmaijer, R. D., & Parr, A. M. 2010, MNRAS, 405, 2439

Packet, W. 1981, A&A, 102, 17

Reggiani, M., & Meyer, M. R. 2013, A&A, 553, A124

Richardson, W. H. 1972, J. Opt. Soc. Am., 62, 55

Rieutord, M. 2016, in Lect. Notes Phys. (Springer Int. Publ.), 914, 101

Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&ARv, 21, 69

Ryans, R. S. I., Dufton, P. L., Rolleston, W. R. J., et al. 2002, MNRAS, 336, 577

Sackmann, I. J. 1970, A&A, 8, 76

Saio, H. 2013, in Lect. Not. Phys. 865, eds. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green (Berlin: Springer Verlag), 159

Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444

Sana, H., Le Bouquin, J.-B., Lacour, S., et al. 2014, ApJS, 215, 15

Semaan, T., Hubert, A. M., Zorec, J., et al. 2013, A&A, 551, A130

Shatsky, N., & Tokovinin, A. 2002, A&A, 382, 92

Shimazaki, H., & Shinomoto, S. 2007, Neural Computation, 19, 1503

Simón-Díaz, S. 2015, in IAU Symp., 307, eds. G. Meynet, C. Georgy, J. Groh, & P. Stee, 194

Simón-Díaz, S., & Herrero, A. 2007, A&A, 468, 1063

Simón-Díaz, S., & Herrero, A. 2014, A&A, 562, A135

Smith, R. C., & Worley, R. 1974, MNRAS, 167, 199

Snow, Jr., T. P. 1981, ApJ, 251, 139

Stoeckley, T. R. 1968, MNRAS, 140, 141

Struve, O. 1931, ApJ, 73, 94

Sundqvist, J. O., Simón-Díaz, S., Puls, J., & Markova, N. 2013, A&A, 559, L10

Tassoul, J.-L. 1978, Theory of rotating stars (Princeton: University Press)

Torres, G. 2010, AJ, 140, 1158

Townsend, R. H. D., Owocki, S. P., & Howarth, I. D. 2004, MNRAS, 350, 189

Tuominen, I. V. 1972, Ann. Acad. Sci. Fennicæ, 392

Unsöld, A., & Struve, O. 1949, ApJ, 110, 455

van Leeuwen, F. 2007, A&A, 474, 653

von Zeipel, H. 1924, MNRAS, 84, 665

Wegner, W. 2006, MNRAS, 371, 185

Wesselius, P. R., van Duinen, R. J., de Jonge, A. R. W., et al. 1982, A&AS, 49, 427

Yudin, R. V. 2001, A&A, 368, 912

Zorec, J. 1986, Ph.D. Thesis, Université Paris VII, Paris

Zorec, J., & Briot, D. 1991, A&A, 245, 150

Zorec, J., & Briot, D. 1997, A&A, 318, 443

Zorec, J., & Royer, F. 2012, A&A, 537, A120

Zorec, J., Divan, L., & Briot, D. 1983, A&A, 126, 205

Zorec, J., Mochkovitch, R., & Divan, L. 1988, C. R. Acad. Sci. Paris, Ser. Sci. Phys., 306, 1265

Zorec, J., Mochkovitch, R. A., & Garcia, A. 1990, in Angular Momentum and Mass Loss for Hot Stars, eds. L. A. Willson, & R. Stalio, NATO ASIC Proc. 316, 239

Zorec, J., Domiciano de Souza, A., & Frémat, Y. 2004, in Stellar Rotation, eds. A. Maeder & P. Eenens, IAU Symp., 215, 21

Zorec, J., Frémat, Y., & Cidale, L. 2005, A&A, 441, 235

Zorec, J., Frémat, Y., Martayan, C., Cidale, L. S., & Torres, A. F. 2007, in Active OB-Stars: Laboratories for Stellar and Circumstellar Physics, eds. A. T. Okazaki, S. P. Owocki, & S. Stefl, ASP Conf. Ser., 361, 539

Zorec, J., Cidale, L., Arias, M. L., et al. 2009, A&A, 501, 297

Zorec, J., Frémat, Y., Domiciano de Souza, A., et al. 2011, A&A, 526, A87

Zorec, J., Frémat, Y., Domiciano de Souza, A., et al. 2016, A&A, submitted (Paper II)
<figure><img src="image_18.png" /><figcaption>**Fig. A.1.** a) Effective temperatures for MK luminosity classes V, IV, and III determined by de Jager & Nieuwenhuijzen (1987) against those extracted from the BCD calibration; b) effective temperatures determined from the uvby photometry using the calibration by Napiwotzki et al. (1993) against those obtained in the BCD system; and c) average effective temperatures determined by Adelman et al. (2002), Nieva (2013), and Fitzpatrick & Massa (2005) against those obtained in the BCD system for a number of common OB stars without emission. The red dashed lines correspond to *x* = *y* diagonal.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. A.2.** Comparison of bolometric luminosities derived using different calibrations of $M_V$ with those derived for the same discrete MK spectral types and luminosity classes using integrated fluxes that are identified here with $\log L(\lambda_1, D)/L_\odot$ (red lines). The identifications are dJM (de Jager & Nieuwenhuijzen 1987), DDS (Deutschman et al. 1976), SK (Aller et al. 1982), and W (Wegner 2006). The blocks **a**), **b**), and **c**) are for average MK luminosity classes V, IV, and III, respectively.</figcaption></figure>

# Appendix A: Fundamental parameters

## A.1. The effective temperature

There are many independent calibrations of effective temperatures given in the literature, which are useful to estimate stellar radii and masses and could be used as references to inquire whether our apparent $T_{\text{eff}}(\lambda_1, D)$ (Zorec et al. 2009) and the critical velocities $V_c$ related to it are affected by systematic deviations. We have chosen: (1) de Jager & Nieuwenhuijzen (1987) because it was used by Cranmer (2005) to discuss the distribution of rotational velocities of Be stars; (2) the calibration of the *ubvy*–Strömgren photometric indices given by Moon & Dworetsky (1985) and revisited by Napiwotzki et al. (1993); and (3) the effective temperatures derived from fitted visible spectral energy distributions by Adelman et al. (2002), Fitzpatrick & Massa (2005), and Nieva (2013). In Fig. A.1a the effective temperatures by de Jager & Nieuwenhuijzen (1987) are compared with the $T_{\text{eff}}(\lambda_1, D)$ for the average luminosity classes V, IV, and III in the stars studied in the present paper. The effective temperatures determined with the *ubvy* indices obtained with the relation given by Napiwotzki et al. (1993) for OBA stars without

emission are compared in Fig. A.1b with $T_{\text{eff}}(\lambda_1, D)$ given in Zorec et al. (2009). Figure A.1c compares the effective temperatures obtained by Adelman et al. (2002) and Nieva (2013) (indicated with AN) with those from Fitzpatrick & Massa (2005) (FM). From these figures we can then conclude that there are no significant systematic deviations between $T_{\text{eff}}(\lambda_1, D)$, $T_{\text{eff}}(\text{AN})$ and $T_{\text{eff}}(\text{FM})$ for stars in the interval of spectral types studied in the present work.

## A.2. The bolometric luminosity

The bolometric luminosities of the program Be stars were determined from integrated monochromatic fluxes. Even though these fluxes were corrected for the presence of circumstellar envelopes/disks, the calculation followed the same rules and data gathered from the same sources as those used in Zorec et al. (2009) for OB stars without emission. Then, we can use the bolometric luminosities for average MK spectral types and luminosity classes that can be determined from this paper to test whether our estimates of bolometric luminosities suffer from some systematic deviations. To this end, we compare our mean
log L<sup>app</sup>/L<sub>☉</sub> with those given by de Jager & Nieuwenhuijzen (1987) and with the bolometric luminosities derived indirectly using calibrations of visual absolute magnitudes M<sub>V</sub>. In this last case, M<sub>V</sub> is transformed into absolute bolometric magnitude using the T<sub>eff</sub>(λ<sub>l</sub>, D) and bolometric corrections (Flower 1996; Torres 2010; Nieva 2013) with M<sub>bol</sub><sup>☉</sup> = +4.742 mag. The calibrations of M<sub>V</sub> against the spectral classes used are from Wegner (2006), Deutschman et al. (1976), and Aller et al. (1982). The log L/L<sub>☉</sub> against the average MK spectral type determined for the average luminosity classes V, IV, and III with the cited methods are shown in Figs. A.2a-c. We can then conclude that our bolometric luminosity estimates are consistent with those determined by other methods. In Fig. A.2, however, there can be a systematic deviation between the BCD luminosities of the average spectral-luminosity classes with those of other authors that stem from the BCD definition of average spectral-luminosity classes. The average log L<sup>app</sup>/L<sub>☉</sub> of the BCD system correspond to the barycenter of the curvilinear quadrilaterals in the (λ<sub>l</sub>, D)-diagram assigned to the average MK spectral types and luminosity classes, while in other calibrations these averages are defined by the specific stellar sample used to carry out the calibrations, which do not necessarily characterize the same centroids as in the BCD system.

## Appendix B: The energy emitted by the star-CE system

The energy emitted by the star-CE system mentioned in Sect. 2.2.2 has been written as (Moujtahid 1998; Moujtahid et al. 2000a)

$$1 + \epsilon_{\lambda} = \alpha_{\lambda}(r, \tau_{\lambda}) + r^{-2}\beta_{\lambda}(r, \tau_{\lambda})[S_{\lambda}(T_{\text{env}})/F_{\lambda}^{*}], \quad (B.1)$$

where $r = R_{*}/R_{\text{env}}$; $R_{\text{env}}$ is an equivalent mean extent of the circumstellar region analyzed, $R_{*}$ is the stellar radius, $F_{\lambda}^{*}$ is the stellar flux, $\tau_{\lambda}$ is the total continuum opacity of the region, and $S_{\lambda}$ is its source function for continuum radiation given here simply by the Planck function in terms of the CE temperature $T_{\text{env}}$; $S_{\lambda} = B_{\lambda}(T_{\text{env}})$. The functions $\alpha_{\lambda}$ and $\beta_{\lambda}$ are given by

$$\alpha_{\lambda}(r, \tau_{\lambda}) = r^{-2} \left[ E_{3}(\tau_{\lambda}) + \mu_{o}^{2} E_{3} \left( \frac{\tau_{\lambda}}{\mu} \right) \right] \\ \beta_{\lambda}(r, \tau_{\lambda}) = 1 - r^{-2} \alpha(r, \tau_{\lambda}) - 2 \mu^{2} E_{3} \left( \frac{2\tau_{\lambda}}{\mu} \right). \quad (B.2)$$

where $E_{3}$ is the exponential integral or order 3 and $\mu_{o} = [1 - r^{2}]^{1/2}$. In the limits for $r \to 0$ and $r \to 1$, the expression in Eq. (B.1) recovers the forms commonly used in the literature

$$1 + \epsilon_{\lambda} = e^{-\tau_{\lambda}} + r^{-2} S_{\lambda} [1 - 2E_{3}(2\tau_{\lambda})] \quad \text{for } r \to 0 \\ 1 + \epsilon_{\lambda} = 2E_{3}(\tau_{\lambda})F_{\lambda} + r^{-2} S_{\lambda} [1 - 2E_{3}(\tau_{\lambda})] \quad \text{for } r \to 1. \quad (B.3)$$

## Appendix C: The equatorial critical linear velocity

In principle, the radiation pressure in the hottest Be stars may produce a significant compensation of the stellar-inside directed effective gravity, so that two critical rotational velocities can be defined (Maeder 2009). Writing the total surface gravity as

$$g_{\text{tot}} = g_{\text{eff}}[1 - \Gamma(\Omega, \theta)], \quad (C.1)$$

<table><caption>**Table B.1.** CE parameters derived for 10 Be stars studied here with the stronger emission in the Hα line.</caption><thead><tr><th>HD</th><th>r ± σ<sub>r</sub></th><th>T<sub>env</sub> ± σ<sub>T</sub></th><th>τ<sub>e</sub> ± σ<sub>e</sub></th><th>τ<sub>V</sub> ± σ<sub>V</sub></th></tr></thead><tbody><tr><td>28497</td><td>0.55 ± 0.04</td><td>13 660 ± 740</td><td>0.79 ± 0.09</td><td>0.06 ± 0.01</td></tr><tr><td>50013</td><td>0.30 ± 0.01</td><td>11 120 ± 100</td><td>0.66 ± 0.03</td><td>0.05 ± 0.01</td></tr><tr><td>68980</td><td>0.36 ± 0.02</td><td>11 060 ± 290</td><td>0.49 ± 0.04</td><td>0.04 ± 0.02</td></tr><tr><td>91465</td><td>0.76 ± 0.27</td><td>10 890 ± 2050</td><td>0.56 ± 0.30</td><td>0.05 ± 0.03</td></tr><tr><td>105435</td><td>0.63 ± 0.01</td><td>14 470 ± 100</td><td>1.42 ± 0.05</td><td>0.05 ± 0.01</td></tr><tr><td>110432</td><td>0.31 ± 0.03</td><td>11 430 ± 990</td><td>0.77 ± 0.14</td><td>0.05 ± 0.01</td></tr><tr><td>120991</td><td>0.83 ± 0.09</td><td>14 280 ± 330</td><td>1.18 ± 0.12</td><td>0.09 ± 0.01</td></tr><tr><td>124367</td><td>0.93 ± 0.04</td><td>15 200 ± 1560</td><td>0.36 ± 0.06</td><td>0.02 ± 0.01</td></tr><tr><td>148184</td><td>0.32 ± 0.01</td><td>11 920 ± 100</td><td>0.91 ± 0.03</td><td>0.06 ± 0.01</td></tr><tr><td>209409</td><td>0.93 ± 0.05</td><td>10 250 ± 340</td><td>0.33 ± 0.14</td><td>0.04 ± 0.02</td></tr></tbody></table>

for which the following definitions hold:

$$g_{\text{eff}} = \frac{GM}{R(\theta)} \{[1 - \eta(\theta) \sin^2 \theta]^2 + \eta^2(\theta) \sin^2 \theta \cos^2 \theta\},$$

$$\eta(\theta) = \frac{\Omega^2 R^3(\theta)}{GM},$$

$$\Gamma(\Omega, \theta) = \frac{\kappa(\theta)L(\Omega)}{4\pi c GM[1 - m(\Omega)]},$$

$$m(\Omega) = \left( \frac{\Omega^2}{2\pi G \bar{\rho}} \right), \quad (C.2)$$

where $g_{\text{eff}}$ is the colatitude-angle dependent surface effective gravity; $R(\theta)$ is the stellar surface radius vector; $\kappa(\theta)$ is the Rosseland mean total opacity for the atmospheric optical depth $\tau = 2/3$; $L(\Omega)$ is the bolometric luminosity changed by the rotational effects, and $\bar{\rho}$ is the stellar density averaged over its entire volume. For $\Gamma(\Omega, \theta) < 1$, the critical velocity is determined by $g_{\text{eff}} = 0$, which for a surface rigid rotation occurs at the equator ($\theta = \pi/2$) and corresponds to $\eta(\theta = \pi/2) = 1.0$ that translates into the well-known expression for the equatorial critical linear velocity

$$V_c = 436.7 \left[ \frac{M/M_{\odot}}{R_c(M, t)/R_{\odot}} \right]^{1/2} \text{ km s}^{-1}. \quad (C.3)$$

In very massive Be stars, i.e., Oe stars with $M \le 40 M_{\odot}$, we could have $\Gamma(\Omega, \theta) = 1$ at some colatitude $\theta$ where $\eta(\theta = \pi/2) < 1.0$. Nevertheless, using the von Zeipel relation to describe the latitudinal variation of the effective temperature, the absorption coefficients calculated for the corresponding gravity and temperature (Kurucz 1993), and knowing that $m(\Omega) \le 0.361$ (Maeder 2009), the condition $\Gamma(\Omega, \theta) = 1$ can be translated into $T_o/T_{\odot} \ge 12.5(g_o/g_{\odot})^{1/4}$ ($T_o$ and $T_{\odot}$ are the effective temperatures of a parent nonrotating star and solar, respectively, while $g_o$ and $g_{\odot}$ are the respective surface gravities). Since of the Be stars in our sample most have $M \le 25 M_{\odot}$ (only 16 of them seem to have $M \ge 25 M_{\odot}$) for which it is $T_o/T_{\odot} \le 12.5(g_o/g_{\odot})^{1/4}$, in our calculation we used the equatorial critical velocity given by Eq. (6).

## Appendix D: Integral equation of probability density distributions

In several sections of this paper we construct integral equations that relate the distributions (frequency of probability densities) of the respective, either true and apparent parameters, or parameters corrected statistically from some physical effect. To
avoid repetitions, we summarize below the basics of the formulation followed in each particular application of the method.

We consider an observed parameter z, whose observed distribution is Z(z) and where z is function of parameters x and y

$$z = z(x, y), \qquad (D.1)$$

with probability densities distributions of X(x) and Y(y), respectively. According to the definition of probability, Z(z) dz must be the surface integral of X(x)Y(y) over the area Δσ between the curves z = constant and z + dz = constant

$$Z(z)dz = \iint_{\Delta\sigma} X(x)Y(y) \, dx \, dy. \qquad (D.2)$$

For some explicit value of x in Eq. (D.2), the width of the strip Δσ is obviously dy = dz/[∂z/∂y] that introduced into Eq. (D.2) leads to

$$Z(z) = \int_{x=x(y)} X(x) \frac{Y[y(x, z)]}{\partial z/\partial y} \, dx, \qquad (D.3)$$

which has to be integrated over the curve x = x(y) given by Eq. (D.1) for the required value of z and within the limits of x established by the specific physical phenomenon analyzed. In the same way, we can obtain and equivalent expression to Eq. (D.3) given by

$$Z(z) = \int_{y=y_c(x)} Y(y) \frac{X[x(y, z)]}{\partial z/\partial x} \, dy. \qquad (D.4)$$

In this paper, we derive the distributions X(x) or Y(y) as counterparts of Z(z) or simply Z(z) once X(x) and Y(y) are defined.

## Appendix E: Determination of pnrc parameters

### E.1. Two-fold iterations

The calculation of parameters $M/M_\odot$, $t/t_{MS}$, $\eta$ and $i$ with Eq. (3) requires two iteration “series” or “stages”. In the first stage, we derive $M/M_\odot$ and $t/t_{MS}$, while in the second stage, we refine the estimates of $\eta$ and $i$.

The models of stellar evolution with rotation used to interpolate masses and ages (Ekström et al. 2012) require the use of fundamental parameters averaged over the rotationally deformed stellar surface. At each n-iteration step of ($M/M_\odot$, $t/t_{MS}$) quantities we establish a list of m-pairs ($\eta_k$, $i_k$), $k = 1, ..., m$ ($m \approx 100$), so that $\eta_{min}(M, t) \le \eta_k \le 1$ and the respective inclinations range in the interval $\pi/2 \ge i_k \ge i_{min}$. From Eq. (3) we draw m-triplets ($T_{\text{eff}}^{\text{pnrc}}(M, t)_k$, $g_{\text{eff}}^{\text{pnrc}}(M, t)_k$, $L_{\text{eff}}^{\text{pnrc}}(M, t)_k$) calculated as $T_{\text{eff}}^{\text{pnrc}}(M, t)_k = T_{\text{eff}}^{\text{app}}/C_T(M, t, \eta_k, i_k)$, $g_{\text{eff}}^{\text{pnrc}}(M, t)_k = g_{\text{eff}}^{\text{app}}/C_G(M, t, \eta_k, i_k)$ and $L_{\text{eff}}^{\text{pnrc}}(M, t)_k = L_{\text{app}}/C_L(M, t, \eta_k, i_k)$, which are then used to determine the m-surface-averaged values $\langle T_{\text{eff}}(M, t, \eta) \rangle_k$, $\langle g(M, t, \eta) \rangle_k$, $\langle L(M, t, \eta) \rangle_k$. Taken by pairs $[\langle T_{\text{eff}}(M, t, \eta) \rangle_k, \langle g(M, t, \eta) \rangle_k]$ and $[\langle T_{\text{eff}}(M, t, \eta) \rangle_k, \langle g(M, t, \eta) \rangle_k]$ as entry parameters to the models of stellar evolution, we determine two independent series of mass and age estimates [$M_k/M_\odot$, $\eta_k$], [$t_k/t_{MS}$, $\eta_k$]. A third series of mass estimates comes from the definition $M_k/M_\odot = (L_k^{\text{pnrc}}/L_\odot)(g_{\text{eff}}^{\text{pnrc}}/g_{\odot,k})(T_\odot/T_{\text{eff},k})^4$. The interpolation of masses and ages is made on the evolutionary models with rotation, labeled with ZAMS rotational velocity previously determined with curves giving the variation of $\Omega/\Omega_c$ with time calculated by Ekström et al. (2012), so that at the

<figure><img src="image_20.png" /><figcaption>Fig. E.1. a) Determination of the pnrc mass $M/M_\odot$ interpolated in stellar evolution models with rotation (black and red lines) and from the relation $M \propto L.g/T_{\text{eff}}^4$ (blue line). b) Determination of $\eta$ and $i$ using the functions $v_{\text{app}}$, $C_T$, $C_L$, and $C_G$ defined in Eq. (3).</figcaption></figure>

iterated stellar age $t_k/t_{MS}$ we recover the $\Omega_k/\Omega_c$ for the given ratio $\eta_k$.

The intersection of the three series of curves [$M_k/M_\odot$, $\eta_k$], [$t/t_{MS}$, $\eta_k$] determine the n-estimate of the stellar mass $M_n/M_\odot$ and age $t_n/t_{MS}$ as shown in Fig. E.1a, which corresponds to the iteration step $n = 1$. These three mass/age curves are only coincident for nonrotating stars. Depending on the star, we need from 5 to 50 iterations to obtain that all three curves in Fig. E.1a determine a unique crossing point.

In the second iteration series, we take advantage of the above estimated $M_n/M_\odot$ and $t_n/t_{MS}$ to establish five relations between $\sin i$ and $\eta$, whose common intersection is considered the refined estimate of $\eta$ and $i$. These curves are the solutions of Eq. (3) written as $\sin i = \sin i(\eta)$ for the given values of $M_n/M_\odot$ and $t_n/t_{MS}$, respectively. A fifth curve is also obtained from the combination $C_R(\eta, i) = C_L^{1/2}/C_T^2$, which may seem redundant, but it is sometimes useful to disentangle the right solution when the intersection of other curves is not unique. Figure E.1b shows the behavior of these relations obtained at the very end of the iteration procedure. This iteration requires another 3 to 15 iteration steps.

The observed parameters [$T_{\text{eff}}^{\text{app}}$, $g_{\text{eff}}^{\text{app}}$, $L_{\text{app}}$, $v_{\text{app}}$] of each program star are considered with their respective uncertainties. So, an observed quantity X is assigned 11 different values within the interval [$X_o - 2.5\sigma_x$, $X_o + 2.5\sigma_x$] using random numbers, which means that for a given star we get 14 641 quartets ($T_{\text{eff}}^{\text{app}}$, $g_{\text{eff}}^{\text{app}}$, $L_{\text{app}}$, $v_{\text{app}}$). Obviously, a certain number of such combinations do not produce any viable solution, however, for most stars we obtain more than $N_{\text{sol}} \approx 10^4$ solutions for [$M/M_\odot$, $t/t_{MS}$, $\eta$, $i$]. In most cases these solutions have asymmetric distributions, which are used as weighting functions to produce weighted averages and their respective equivalent $1\sigma$ dispersions.
<table><caption>**Table E.1.** Test *pnrc* solutions parameters.</caption><thead><tr><th rowspan="2">η</th><th rowspan="2">i°</th><th colspan="4">Input apparent model-data</th><th colspan="7">Solutions for pnrc parameters</th></tr><tr><th>T<sub>eff</sub><br>K</th><th>log g</th><th>log L/L<sub>☉</sub></th><th>Vsin i</th><th>Vsin i<br>corr</th><th>η</th><th>i°</th><th>T<sub>eff</sub><sup>°</sup><br>K</th><th>log g<sub>o</sub></th><th>log L<sub>☉</sub>/L<sub>☉</sub></th><th>Vsin i<br>corr</th><th>V<sub>c</sub></th><th>V<sub>c</sub> km s<sup>-1</sup></th></tr></thead><tbody><tr><td>0.00</td><td>0°</td><td>20 000</td><td>4.00</td><td>3.461</td><td>0</td><td>0</td><td>0.00</td><td>0°</td><td>20 000</td><td>4.00</td><td>3.461</td><td>0</td><td>462</td><td></td></tr><tr><td>0.20</td><td>20°</td><td>19 665</td><td>3.95</td><td>3.496</td><td>81</td><td>81</td><td>0.20</td><td>20°</td><td>20 001</td><td>4.00</td><td>3.461</td><td>81</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.22 ± 0.17</td><td>23° ± 10</td><td>20 093 ± 79</td><td>4.01 ± 0.02</td><td>3.462 ± 0.025</td><td>82 ± 6</td><td>465 ± 5</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.20</td><td>50°</td><td>20 002</td><td>4.00</td><td>3.461</td><td>182</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.20 ± 0.13</td><td>48° ± 11</td><td>19 915 ± 79</td><td>3.99 ± 0.03</td><td>3.444 ± 0.018</td><td>184 ± 14</td><td>463 ± 5</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.20</td><td>80°</td><td>19 999</td><td>4.00</td><td>3.461</td><td>234</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.22 ± 0.09</td><td>60° ± 10</td><td>19 792 ± 338</td><td>3.98 ± 0.03</td><td>3.424 ± 0.039</td><td>237 ± 18</td><td>462 ± 6</td><td></td></tr><tr><td>0.50</td><td>20°</td><td>19 555</td><td>3.90</td><td>3.572</td><td>119</td><td>122</td><td>0.46</td><td>22°</td><td>20 116</td><td>3.99</td><td>3.488</td><td>122</td><td>461</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.45 ± 0.15</td><td>24° ± 6</td><td>20 114 ± 75</td><td>3.98 ± 0.02</td><td>3.496 ± 0.015</td><td>122 ± 9</td><td>460 ± 5</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.45</td><td>53°</td><td>20 092</td><td>4.00</td><td>3.476</td><td>270</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.45 ± 0.09</td><td>54° ± 6</td><td>20 110 ± 224</td><td>4.00 ± 0.01</td><td>3.479 ± 0.017</td><td>271 ± 21</td><td>462 ± 5</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.51</td><td>71°</td><td>19 811</td><td>3.95</td><td>3.426</td><td>355</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.54 ± 0.11</td><td>74° ± 8</td><td>20 051 ± 598</td><td>3.96 ± 0.02</td><td>3.448 ± 0.062</td><td>363 ± 36</td><td>466 ± 11</td><td></td></tr><tr><td>0.75</td><td>20°</td><td>19 143</td><td>3.89</td><td>3.615</td><td>126</td><td>143</td><td>0.75</td><td>20°</td><td>20 004</td><td>4.00</td><td>3.461</td><td>143</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.73 ± 0.13</td><td>20° ± 2</td><td>19 965 ± 108</td><td>4.00 ± 0.00</td><td>3.457 ± 0.012</td><td>141 ± 11</td><td>462 ± 1</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.75</td><td>50°</td><td>20 000</td><td>4.00</td><td>3.460</td><td>320</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.76 ± 0.11</td><td>50° ± 4</td><td>19 967 ± 262</td><td>4.00 ± 0.01</td><td>3.454 ± 0.026</td><td>317 ± 24</td><td>462 ± 4</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.75</td><td>80°</td><td>19 999</td><td>4.00</td><td>3.461</td><td>411</td><td>462</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.76 ± 0.14</td><td>78° ± 7</td><td>19 758 ± 536</td><td>3.99 ± 0.02</td><td>3.434 ± 0.049</td><td>403 ± 42</td><td>459 ± 10</td><td></td></tr></tbody></table>

This procedure enables us to obtain Stoeckley’s correction $\Sigma(M, t, \eta, i)$, which is consistent with the stellar fundamental parameters.

## E.2. Reliability of the inferred pnrc parameters

We tested the procedure described in Sects. 3.5.1 and E.1 using several model stars, where the apparent parameters were simulated with FASTROT using model atmospheres with solar metallicity as required for the program stars. We performed many tests, in particular for a typical model star characterized by $T^{pnrc} = 20 000 \pm 500$ K, $\log g^{pnrc} = 4.0 \pm 0.01$, $\log L^{pnrc}/L_{\odot} = 3.461 \pm 0.025$, so that $M/M_{\odot} = 7.31 \pm 0.08$, $t/t_{MS} = 0.710 \pm 0.016$ and $V^{pnrc} = 462$ km s$^{-1}$. This mathematical object was assumed to rotate successively at $\eta = 0.2, 0.5$ and $0.75$, and each time seen at inclinations $i = 20^{\circ}, 50^{\circ}$ and $80^{\circ}$. The results obtained for this model object are given in Table E.1, where we can see that the algorithm recovers the right pnrc parameters if they are treated without uncertainties. However, when the input parameters are taken with their uncertainties, as occurs with the observed stars, unavoidable deviations appear because the error propagation is asymmetric. Nevertheless, in Table E.1 we see that the error bars of the recovered pnrc quantities are not larger than the input uncertainties quoted above in this paragraph, which warrants that the algorithm is consistent.

## Appendix F: Approximate or auxiliary distributions

Let $x$ be a random variable defined in the interval $0 \le x \le 1$, whose probability density distribution is $\Pi(x)$. When the moments of $\Pi(x)$ are known up to degree $n = 3$, it is possible to derive simplified analytical expressions $\bar{\Pi}(x)$ that can be considered statistically quasi-equivalent to $\Pi(x)$ because they reproduce the same average value of $x$, the same variance (or

<table><caption>**Table F.1.** First four moments of the function $\Pi(x)$ and the parameters of its quasi-equivalent trapezoidal and triangular functions $\bar{\Pi}(x)$.</caption><thead><tr><th>γ</th><th>⟨x⟩</th><th>⟨x<sup>2</sup>⟩</th><th>⟨x<sup>3</sup>⟩</th><th>⟨x<sup>4</sup>⟩</th></tr></thead><tbody><tr><td>-0.75</td><td>0.8313</td><td>1.0000</td><td>1.1687</td><td>1.4643</td></tr><tr><td>0.00</td><td>0.8875</td><td>1.2250</td><td>1.5625</td><td>2.4551</td></tr><tr><td>0.75</td><td>1.0897</td><td>1.6750</td><td>2.2602</td><td>3.4459</td></tr></tbody><thead><tr><th colspan="5">Trapezoidal Π̄(x)</th></tr><tr><th>γ</th><th>a</th><th>b</th><th>p</th><th>q</th><th>⟨x<sup>4</sup>⟩ p%</th></tr></thead><tbody><tr><td>-0.75</td><td>0.0756</td><td>0.8676</td><td>-0.6762</td><td>0.9503</td><td>1.4465 36.3</td></tr><tr><td>0.0</td><td>0.0892</td><td>0.9108</td><td>0.0</td><td>0.6086</td><td>2.4427 66.5</td></tr><tr><td>0.75</td><td>0.1324</td><td>0.9244</td><td>0.6762</td><td>0.2741</td><td>3.4442 36.3</td></tr></tbody><thead><tr><th colspan="5">Triangular Π̄(x)</th></tr><tr><th>γ</th><th>p<sub>1</sub></th><th>q<sub>1</sub></th><th>p<sub>2</sub></th><th>q<sub>2</sub></th><th>⟨x<sup>4</sup>⟩ p%</th></tr></thead><tbody><tr><td>-0.75</td><td>2.5554</td><td>0.1600</td><td>-1.3331</td><td>1.3408</td><td>1.4689 99.7</td></tr><tr><td>0.00</td><td>1.4816</td><td>0.1199</td><td>-1.4816</td><td>1.6015</td><td>2.4714 97.1</td></tr><tr><td>0.75</td><td>1.3331</td><td>0.0077</td><td>-2.5554</td><td>2.7154</td><td>3.4508 99.7</td></tr></tbody></table>

dispersion), and the same skewness (or symmetry) of the $\Pi(x)$. The degree $n$ can be higher than 3, so that moments of higher order can accounted for. We limit our discussion to $n = 3$ because we compare the Lucy-Richardson deconvolution method to the Cranmer method, where $n = 3$ was used. Thus, we develop here simple rules to derive trapezoidal or truncated linear functions, and triangular quasi-equivalent distribution functions $\bar{\Pi}(x)$.

### F.1. Trapezoidal functions

We assume that $\bar{\Pi}(x)$ is given by

$$ \bar{\Pi}(x) = \begin{cases} px + q & x \in [a, b] \\ 0 & x \notin [a, b], \end{cases} \qquad (F.1) $$
where $p$ and $q$, as well as the boundaries $a$ and $b$ are parameters that we need to determine. We then ask that the moments are conserved, i.e.,

$$ \langle x^n \rangle = \int_0^1 x^n \Pi(x) \, dx = \int_a^b x^n \bar{\Pi}(x) \, dx, \qquad (F.2) $$

where $n = 0$ indicates the normalization condition for both $\Pi(x)$ and $\bar{\Pi}(x)$. Calling $\zeta = a + b$ and $\xi = a - b$, the normalization of $\bar{\Pi}(x)$ and its moment of first order requires that

$$ \begin{aligned} \frac{1}{2} p \zeta + q &= \frac{1}{\xi} \\ \frac{1}{6} p \xi^3 + \zeta &= 2 \langle x \rangle, \end{aligned} \qquad (F.3) $$

while the moments of order $n = 2$ and $3$ lead to

$$ \begin{aligned} \lambda^3 + 3P\lambda + 2Q &= 0 \\ \zeta^2 - 4\langle x \rangle\zeta + 4\langle x^2 \rangle &= \frac{\xi^2}{3}, \end{aligned} \qquad (F.4) $$

where the following definitions hold:

$$ \begin{aligned} \lambda &= \zeta - 2\langle x \rangle \\ P &= 2\langle x \rangle^2 - 2\langle x^2 \rangle \\ Q &= 10\langle x \rangle^3 - 15\langle x^2 \rangle\langle x \rangle + 5\langle x^3 \rangle. \end{aligned} \qquad (F.5) $$

Solving in Eq. (F.4) the algebraic equation of 3rd order for $\lambda$, we derive an expression for $\zeta$, which in turn determines $\xi$ using the 2nd relation in Eq. (F.4). The constants $a$ and $b$ are then obtained from Eq. (F.3), which lead to $p$ and $q$.

## F.2. Triangular functions

We assume that $\bar{\Pi}(x)$ is defined as

$$ \bar{\Pi}(x) = \begin{cases} \bar{\Pi}_1(x) = p_1 x + q_1 & x \in [a, c] \\ \bar{\Pi}_2(x) = p_2 x + q_2 & x \in [c, b], \end{cases} $$

where the constants to be determined are $p_1$, $p_2$, $q_1$, and $q_2$. The conditions $\bar{\Pi}_1(a) = \bar{\Pi}_2(b)$ and $\bar{\Pi}_1(c) = \bar{\Pi}_2(c)$ produce

$$ \begin{aligned} p_1 a + q_1 &= 0 \\ p_2 b + q_2 &= 0 \\ (p_1 - p_2)c &= q_2 - q_1. \end{aligned} \qquad (F.6) $$

The normalization of $\bar{\Pi}(x)$ imposes

$$ \begin{aligned} p_1(b-a)(c-a) &= 2 \\ p_2(b-a)(c-b) &= 2. \end{aligned} \qquad (F.7) $$

We then impose that the moments $\langle x^n \rangle$ ($n = 1$ to 3) are conserved as in Eq. (F.2), so that

$$ \begin{aligned} A &= 3\langle x \rangle = a + b + c \\ B &= 9\langle x \rangle^2 - 6\langle x^2 \rangle = ab + ac + bc \\ C &= 10\langle x^3 \rangle + 27\langle x \rangle^3 - 36\langle x \rangle\langle x^2 \rangle = abc, \end{aligned} $$

which lead to the following 3rd order equation that enables us to calculate the boundary $a$

$$ a^3 - a^2 A + aB - C = 0. \qquad (F.8) $$

Defining $\alpha = b + c = A - a$ and $\beta = bc = C/a$, we obtain the 2nd order algebraic equation for the boundary $b$,

$$ b^2 - \alpha b + \beta = 0. \qquad (F.9) $$

<figure><img src="image_21.png" /><figcaption>Fig. F.1. Left side: distribution functions $\Pi(x)$ (thick black curves) and their quasi-equivalent trapezoidal (blue curves) and triangular (red curves) $\bar{\Pi}(x)$ distributions up to the third moment. Right side: respective cumulative distributions $\pi(x) = \int_0^x \Pi(z) \, dz$.</figcaption></figure>

Once $a$ and $b$ are known, we get $c = \alpha - b$. Then, from Eq. (F.7) we readily obtain $p_1$ and $p_2$, so that $q_1$ and $q_2$ can be drawn from Eq. (F.6).

For the sake of an exercise, we adopt the following distribution function$^3$:

$$ \Pi(x) = \left\{ \begin{aligned} &0.3183 [1 - (2x - 1)^2]^{1/2} \\ &+ 0.375 [1 - (2x - 1)^2] \end{aligned} \right\} \times [1 + \gamma(2x - 1)], $$

defined for $0 \le x \le 1$, which is shown in Fig. F.1, for $\gamma = -0.75$, 0.0 and 0.75. From the relations obtained in Sects. F.1 and F.2 we derive the parameters of the quasi-equivalent trapezoidal and triangular distributions $\bar{\Pi}(x)$ given in Table F.1, and the respective functions $\bar{\Pi}(x)$ are shown in Fig. F.1 (blue and red curves). In this figure the respective cumulative distributions $\pi(x) = \int_0^x \Pi(z) \, dz$ (right side blocks) are also shown.

For the approximations $\bar{\Pi}(x)$ obtained here, where by definition only the first three moments are conserved, the most important question to solve is, which of these moments should be considered the closer parent distribution to the actual $\Pi(x)$ one? To this end we added in Table F.1 the moments of order four $\langle x^4 \rangle$ issued from the respective simplified functions $\bar{\Pi}(x)$ and the probabilities $p$% that, according to the Kolmogorov-Smirnov test, can be considered parent with $\Pi(x)$. We thus immediately note that in Cranmer's attempts the triangular would clearly be a better choice than the trapezoidal. Nevertheless, depending on

<footnote>$^3$ This function is inspired on the classical broadening function of spectral lines by rotation. There is no specific reason to use this expression to do the exercise.</footnote>
the statistical problem, the suited simplified functional form $\bar{\Pi}(x)$
is probably different. Moreover, this leads us to conclude that it
is not enough to claim that Gaussian or Maxwellian distributions
are consistent with the observed distributions simply because the
first two moments (average and dispersion) coincide, as the lit-
erature sometimes concludes (e.g., Mokiem et al. 2006; Dufton
et al. 2006b; Hunter et al. 2008).

The approximate distributions are useful when Stoeckley’s
correction can be neglected. Nevertheless, using the relations

between moments of apparent and true distributions, $\Psi(v)$ and
$\Phi(u)$ given by Chandrasekhar & Münch (1950), provide a first
good guess for the searched parameters defining the functions
$\bar{\Phi}(u)$ that are supposed to represent $\Phi(u)$. Trapezoidal and trian-
gular distributions $\bar{\Phi}(u)$ are dependent on four independent pa-
rameters, which means that apart from the normalization con-
dition, they can reproduce only the first three moments of the
actual distribution $\Phi(u)$.