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

# Properties of Galactic B[e] Supergiants. V. 3 Pup–Constraining the Orbital Parameters and Modeling the Circumstellar Environments

A. S. Miroshnichenko<sup>1,2,3</sup>①, S. Danford<sup>1</sup>, S. V. Zharikov<sup>4</sup>, V. G. Klochkova<sup>5</sup>, E. L. Chentsov<sup>5</sup>, D. Vanbeveren<sup>6</sup>, O. V. Zakhozhay<sup>7</sup>, N. Manset<sup>8</sup>, M. A. Pogodin<sup>2</sup>, C. T. Omarov<sup>3,9</sup>, A. K. Kuratova<sup>10</sup>, and S. A. Khokhlov<sup>10</sup>②

<sup>1</sup>Department of Physics and Astronomy, University of North Carolina at Greensboro, P.O. Box 26170, Greensboro, NC 27402-6170, USA

<sup>2</sup>Pulkovo Astronomical Observatory of the Russian Academy of Sciences, Pulkovskoe shosse 65-1, St. Petersburg, 196140, Russia

<sup>3</sup>Fesenkov Astrophysical Institute, Observatory, 23, Almaty, 050020, Kazakhstan

<sup>4</sup>Instituto de Astronomia, Universidad Nacional Autónoma de México, Apdo. Postal 877, Ensenada, 22800, Baja California, México

<sup>5</sup>Special Astrophysical Observatory of the Russian Academy of Sciences, Nizhnyj Arkhyz, Zelenchukskiy region, Karachai-Cherkessian Republic, 369167, Russia

<sup>6</sup>Astronomy and Astrophysics Research Group, Vrije Universiteit Brussel, Pleinlaan 2, B-1050 Brussels, Belgium

<sup>7</sup>Main Astronomical Observatory, National Academy of Sciences of Ukraine, Kyiv 03680, Ukraine

<sup>8</sup>CFHT Corporation, Kamuela, HI 96743, USA

<sup>9</sup>National Center of Space Research and Technology, Shevchenko St. 15, Almaty, 050010, Kazakhstan

<sup>10</sup>Al-Farabi Kazakh National University, Al-Farabi Ave. 71, 050040, Almaty, Kazakhstan

Received 2020 March 17; revised 2020 May 4; accepted 2020 May 15; published 2020 July 1

## Abstract

We report the results of long-term spectroscopic monitoring of the A-type supergiant with the B[e] phenomenon 3 Pup = HD 62623. We confirm earlier findings that it is a binary system. The orbital parameters were derived using cross-correlation of the spectra in a range of 4460–4632 Å, which contains over 30 absorption lines. The orbit was found to be circular with a period of $137.4 \pm 0.1$ days, radial velocity semiamplitude $K_1 = 5.0 \pm 0.8$ km s$^{-1}$, systemic radial velocity $\gamma = +26.4 \pm 2.0$ km s$^{-1}$, and mass function $f(m) = (1.81_{-0.76}^{+0.97}) \times 10^{-3} M_\odot$. The object may have evolved from a pair with initial masses of $\sim 6.0 M_\odot$ and $\sim 3.6 M_\odot$ with an initial orbital period of $\sim 5$ days. Based on the fundamental parameters of the A-supergiant (luminosity $\log L/L_\odot = 4.1 \pm 0.1$ and effective temperature $T_{\text{eff}} = 8500 \pm 500$ K) and evolutionary tracks of mass-transferring binaries, we found current masses of the gainer $M_2 = 8.8 \pm 0.5 M_\odot$ and donor $M_1 = 0.75 \pm 0.25 M_\odot$. We also modeled the object’s IR-excess and derived a dust mass of $\sim 5 \times 10^{-5} M_\odot$ in the optically thin dusty disk. The orbital parameters and properties of the Hα line profile suggest that the circumstellar gaseous disk is predominantly circumbinary. The relatively low mass of the gainer led us to a suggestion that 3 Pup should be excluded from the B[e] supergiant group and moved to the FS CMa group. Overall these results further support our original suggestion that FS CMa objects are binary systems, where an earlier mass transfer caused formation of the circumstellar envelope.

Unified Astronomy Thesaurus concepts: Emission line stars (460); Binary stars (154); High resolution spectroscopy (2096)

Supporting material: machine-readable table

## 1. Introduction

The object 3 Pup (HD 62623, HR 2996) is the brightest ($V \sim 4.0$ mag) among those exhibiting the B[e] phenomenon. The latter is defined as the presence of emission lines, including forbidden, in the spectra of B-type stars as well as an IR excess due to dust radiation. The emission-line spectrum of the star was first detected by Merrill (1934). The first identification of spectral lines in a wavelength range from 3300 Å to $H\alpha$ was published by Swings (1950). A more detailed history of spectral studies of 3 Pup along with a line list found at a high resolution ($R = 60,000–70,000$) in a range from 3682 to 8863 Å can be found in Chentsov et al. (2010). An atlas of the spectrum of 3 Pup in the spectral range from 3920 Å to 6920 Å at $R = 60,000$ is presented by Klochkova et al. (2015). Parts of the spectrum at longer wavelengths from 7280 Å to 7340 Å and 8480 Å to 8680 Å at $R = 15,000–18,000$ are presented by Aret et al. (2016).

Regular radial velocity (RV) variations of the absorption lines in the spectrum of 3 Pup were first detected by Johnson & Neubauer (1946), who reported an orbital period of 137.767 days and an RV semiamplitude of $3.60 \pm 0.45$ km s$^{-1}$. These photographic data with the addition of several higher-resolution CCD spectra were reanalyzed by Plets et al. (1995). Using a

Fourier analysis technique, these authors found two equally probable orbital periods of 138.5 and 161.1 days. They also estimated the object’s distance modulus to be 10.8–12.4 mag corresponding to distances from 1.4 to 3.0 kpc. This result led Plets et al. (1995) to a very large luminosity of the star ranging from $(1.1 \pm 0.5) \times 10^5 L_\odot$ and a conclusion that 3 Pup was a massive star. Such a large distance and a high luminosity were not supported by the results of Chentsov et al. (2010), who derived a distance of 650 pc from the spectral classification of the star (A2.7 lb) and comparison of its spectrum with those of other early A-type supergiants.

The IR excess of 3 Pup was analyzed by Rovero & Ringuelet (1994), Plets et al. (1995), Stee et al. (2004), and Meilland et al. (2010). The former authors attributed the excess to chromospheric emission, while the others clearly showed that it is due to radiation of the circumstellar dust. A strong emission band, indicative of optically thin radiation of silicate dusty particles was detected in the low-resolution IRAS spectrum. Stee et al. (2004) calculated a model SED of the star and the gas-and-dust disk-like envelope using a stellar temperature of 10,000 K, a typical accelerated stellar wind of a single supergiant with a latitudinal mass flux distribution, and spherical silicate dusty grains. The best model found by these authors contained dusty
<table><caption>**Table 1**
Observatories and Instruments Used for Spectroscopic Observations</caption><thead><tr><th>Observatory ID</th><th>Telescope</th><th>Instrument</th><th>Resolution λ/Δλ</th><th>Location</th></tr></thead><tbody><tr><td>TCO</td><td>0.81 m at the Three College Observatory</td><td>eShel (Shelyak Instruments)<sup>a</sup></td><td>12,000</td><td>North Carolina, USA</td></tr><tr><td>SAO</td><td>6 m at the Special Astrophysical Observatory of the Russian Academy of Sciences</td><td>NES (Panchuk et al. 2017)</td><td>60,000</td><td>Nizhniy Arkhyz, Russia</td></tr><tr><td>OAN SPM</td><td>2.1 m at the Observatorio Astronómico Nacional San Pedro Martir</td><td>REOSC</td><td>18,000</td><td>Baja California, Mexico</td></tr><tr><td>McD</td><td>2.7 m Harlan J. Smith at the McDonald Observatory</td><td>Tull coudé TS2 (Tull et al. 1995)</td><td>60,000</td><td>Mt. Locke, Texas, USA</td></tr><tr><td>CFHT</td><td>3.6 m Canada-France-Hawaii Telescope</td><td>ESPaDOnS (Manset & Donati 2003)</td><td>65,000</td><td>Maunakea, Hawaii, USA</td></tr></tbody></table>

**Note.**

grains with a diameter of 1 µm and was consistent with near-IR interferometric data.

More recent studies (e.g., Monnier et al. 2004; Meilland et al. 2010; Millour et al. 2011) reported IR interferometry of the object, resolved the dusty disk around it, and concluded that the material accumulation in its equatorial plane is most likely due to the presence of a low-mass component. Also, spectroscopic studies by Kraus et al. (2015) and Aret et al. (2016) in the visual and near-IR regions suggested a Keplerian motion of the material in the circumstellar disk around the A-supergiant.

Nevertheless, there is still much uncertainty in the stellar and circumstellar parameters of the object. In particular, there is no generally accepted agreement on the luminosity of the A-supergiant, mass function of the binary, and even the orbital period. In this paper we report the results of a long-term high-resolution spectroscopic monitoring of 3 Pup and a new attempt to model the dusty portion of the system IR excess with the goal to further constrain the system parameters.

## 2. Observations

Spectroscopic observations of 3 Pup were obtained by A.M. and S.D. in 2012–2020 at the Three College Observatory in North Carolina (131 spectra). Additionally 17 spectra were obtained in 2004–2018 at four additional observatories. Table 1 summarizes these observations.

The data obtained at OAN SPM, McDonald (McD), and TCO were reduced in a standard way using the *echelle* task of IRAF. Observations obtained at CFHT were reduced with the Upena and Libre-ESPRIT software packages (Donati et al. 1997). Typical uncertainties in the wavelength calibration are <0.5 km s⁻¹ for the CFHT, McDonald, and TCO data and ~1 km s⁻¹ for OAN SPM. RV standard stars were observed every night at TCO and SPM to control the wavelength calibration. Unexpected deviations from the regular RV variations described in Section 3.1 were detected only several times and might have been due to either the pulsational activity of the A-supergiant or to effects of flexure (the spectrograph at OAN SPM is mounted at the Cassegrain focus). Log of the observations along with the measurements of the Hα line properties and cross-correlated RV is presented in Table 2.

## 3. Data Analysis

### 3.1. Binary Orbit

The optical spectrum of 3 Pup and line identification were described in detail in Chentsov et al. (2010). Our goal in this study was to investigate long-term line position and profile variations using a much larger data set than those reported previously. To reach this goal, we have measured the emission peaks intensity ratios in the Hα profile and positions of over 30 absorption lines of Fe II, Mg II, Ti II, and Cr II in a blue region from 4460 Å to 4632 Å as well as in a red region around the Si II lines at 6347 Å and 6371 Å. Spectra with both regions of absorption lines were normalized to a local continuum and cross-correlated against a template spectrum, which was chosen to be a TCO spectrum of 3 Pup obtained on 2017 February 17 (signal-to-noise ratio in continuum over 200, see Figure 1), using the *rvsao* package in IRAF. The template spectrum was obtained near an orbital phase, at which the observed RV is close to the systemic velocity. Also, using an individual spectrum as a template is better than an averaged one to avoid a spectral line broadening because of additional variations due to other processes, such as pulsations (see a brief discussion below).

The spectrum of HD 55036 shown for comparison in Figure 1 was selected from a sample of A2–A3 Ib supergiants that have been observed at TCO and OAN SPM. We show it to demonstrate that the spectrum of 3 Pup looks like that of a typical supergiant in the regions with no line emission.

The main result found from the cross-correlation is a very well-defined pattern of the RV variations (see Figure 2). It turned out to be sinusoidal within the measurement uncertainties, indicating a circular orbit of the A-supergiant around the center of mass of the system. The best fit to the blue region data has the following parameters: the orbital period $P_{orb} = 137.4 \pm 0.1$ days, the RV semiamplitude $K = 5.0 \pm 0.8$ km s⁻¹, and the epoch of the superior conjunction $HJD_0 = HJD245,3325.7 \pm 3.0$. The systemic velocity $\gamma = +26.4 \pm 2.0$ km s⁻¹ was determined by adding the average RV from cross-correlation and the average RV of all the absorption lines in the blue region of the template spectrum.

<footnote>a http://www.shelyak.com</footnote>
<table><caption>**Table 2**
Log of Spectroscopic Observations of 3 Pup</caption><thead><tr><th>Date</th><th>HJD-2450,000</th><th>Obs.</th><th>Range, Å</th><th>Phase</th><th>RV</th><th>σ(RV)</th><th>I<sub>V</sub></th><th>I<sub>R</sub></th><th>I<sub>d</sub></th><th>V/R</th></tr></thead><tbody><tr><td>2004 Dec 25</td><td>3365.042</td><td>CFHT</td><td>3600–10500</td><td>0.287</td><td>-5.5</td><td>0.3</td><td>1.92</td><td>2.84</td><td>0.84</td><td>0.68</td></tr><tr><td>2006 Dec 12 (b)</td><td>4082.943</td><td>SPM</td><td>3850–6865</td><td>0.511</td><td>...</td><td>...</td><td>1.95</td><td>2.31</td><td>1.24</td><td>0.84</td></tr><tr><td>2007 Nov 15</td><td>4420.016</td><td>SPM</td><td>3600–6775</td><td>0.964</td><td>-1.1</td><td>0.4</td><td>1.63</td><td>2.12</td><td>1.36</td><td>0.77</td></tr><tr><td>2007 Nov 20 (b)</td><td>4424.982</td><td>SPM</td><td>3600–6775</td><td>0.000</td><td>...</td><td>...</td><td>1.65</td><td>2.30</td><td>1.48</td><td>0.72</td></tr><tr><td>2008 Oct 10</td><td>4750.017</td><td>SPM</td><td>3690–6800</td><td>0.365</td><td>-4.4</td><td>0.5</td><td>1.75</td><td>2.98</td><td>1.15</td><td>0.59</td></tr><tr><td>2008 Nov 4 (a)</td><td>4774.627</td><td>SAO</td><td>4462–5926</td><td>0.545</td><td>-0.5</td><td>0.1</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>2008 Dec 11 (c)</td><td>4811.881</td><td>McD</td><td>3600–10140</td><td>0.816</td><td>...</td><td>...</td><td>1.60</td><td>2.88</td><td>1.14</td><td>0.56</td></tr></tbody></table>

**Note.** Log of spectroscopic observations of 3 Pup. The full table is shown in the electronic version of the paper. Column information: (1) calendar date (YY MM DD), (2) Julian date (JD-2450,000), (3) observatory ID (see Table 1), (4) spectral range observed, (5) orbital phase according to the RV solution (see the text), (6) radial velocity in kilometers per second derived by cross-correlation in the range 4460–4632 Å (see Figure 1), (7) uncertainty in the radial velocity determination in kilometers per second, and (8)–(11) parameters of the Hα line profiles—blue peak intensity in continuum units ($I_V$), red peak intensity ($I_R$), intensity of the central depression ($I_d$), and the peak intensity ratio ($V/R$). Comments on the spectra with no RV measurements: (a) a large portion or the entire cross-correlation region was not observed; (b) a systematic error in the wavelength calibration; and (c) the region damaged by reflection on the CCD chip.

(This table is available in its entirety in machine-readable form.)

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Template for cross-correlation of the absorption spectrum of 3 Pup taken at TCO (lower spectrum). Part of the spectrum of HD 55036, an A3 Ib supergiant, taken at the 2 m OAN SPM telescope and convolved with the TCO spectral resolution is shown for comparison. Intensity is normalized to the local continuum, the wavelength scale is shifted to match the line positions and given in angstroms. Identifications of some lines are shown. Most other lines belong to Fe II.</figcaption></figure>

Although the RV measurement accuracy is of the order of 0.3 km s⁻¹, the uncertainty of the derived RV semiamplitude due to the orbital motion is nearly three times higher. It is most likely due to short-term variations of the line positions mentioned by Chentsov et al. (2010). Such variations may be caused by pulsations typical of early-type supergiants (e.g., Saio et al. 2013).

The RV curve folded with the orbital period is shown in the lower panel of Figure 2. The cross-correlation of the spectra in the red region showed virtually the same results as for the blue region. The orbital period and the RV semiamplitude result in a mass function $f(m) = (1.81^{+0.97}_{-0.76}) \times 10^{-3} M_\odot$, which is discussed in Section 4.

## 3.2. Other Spectral Variations

Previous papers that reported high-resolution spectra of 3 Pup contained only limited discussions of the line profile variations because of the small number of spectra obtained. For example, Chentsov et al. (2010) mentioned profile variations of

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

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Upper panel: Fourier power spectrum for the RV variations. Lower panel: RV phase curve based on cross-correlation of the 4460–4632 Å spectral region. The RV scale is heliocentric and shown in kilometers per second. Filled circles represent TCO data, open circles show SAO data, diamonds show SPM data, and the triangle shows CFHT data. Bottom section of the plot shows deviations of the observed RV from the calculated one.</figcaption></figure>

Fe II and Na I lines comparing only two spectra. Our collection allows tracing variations of spectral lines over a period of nearly two decades. Overall, the absorption-line spectrum of 3 Pup in the visual region is typical of Ib supergiants, and
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Hα line profiles in the SAO spectrum taken on 2013 October 14 (solid line) and the McDonald spectrum taken on 2008 December 11 (dashed line). Both spectra have the same R = 60,000. Telluric lines were not removed from the data. Intensity is normalized to the local continuum, the RV is shown in kilometers per second.</figcaption></figure>

absorption lines show very weak variations of their profiles, taking into account uncertainties due to continuum normalization and wavelength calibration.

We analyzed the Hα line double-peaked profile variations and found several features of its long-term behavior. The first one concerns the strength of the central depression, which varied between 1.08 and 1.48 of the continuum intensity even in our highest-resolution spectra in 2006–2012. Since 2013 October it has never exceeded 1.03 of the continuum intensity. As seen in Figure 3, the blueshifted peak was lower during the former period compared to the latter one. Also, the blueshifted peak was accompanied by a weaker and even bluer-shifted peak at an RV of ~-100 km s⁻¹.

Other features show correlation with the orbital phase (see Figure 4). The peak intensity ratio V/R shows a tendency for a lower V/R near inferior conjunction (orbital phase = 0.5), when the A-supergiant is located in front of the secondary component. At the same time, the red emission peak is noticeably stronger around this phase. The values of V/R and the red peak intensity vary from cycle to cycle, but overall the dependence keeps its shape. Also, the emission peak separation and the red peak position (measured by fitting an area around the peak with a Gaussian) exhibit a well-defined maximum near the same phase.

### 3.3. Spectral Energy Distribution

The IR excess in the SED of 3 Pup has been discussed in several papers. Plets et al. (1995) mentioned that the dusty envelope extends from R<sub>in</sub> = 35R<sub>*</sub> to R<sub>out</sub> = 5000R<sub>*</sub>, where R<sub>*</sub> is the radius of the A-supergiant, with a radial density profile of ρ ∝ r<sup>-1.3</sup> but have not given the modeling details. The best model of Stee et al. (2004) has similar parameters (R<sub>in</sub> = 20R<sub>*</sub>, R<sub>out</sub> = 5000R<sub>*</sub>, and ρ ∝ r<sup>-2.0</sup>), using a grain size of 1 µm and IR fluxes at wavelengths up to 100 µm. The latter authors and Meilland et al. (2010) focused on explaining the mid-IR interferometry results.

We tried reproducing the SED shape in a wider wavelength range. The SED of 3 Pup was constructed from various sources, which included UV, optical photometry, near-IR, and IRAS photometry (presented in Plets et al. 1995), fluxes from

more recent IR surveys (Wide-field Infrared Survey Explorer, AKARI), and a set of submillimetric fluxes from Jura et al. (2001). The star has been found to show little variability in the visual spectral range (0.03 mag, e.g., Sterken et al. 1993; Adelman & Albayrak 1997). Therefore, combining data taken at different times to construct the SED is justified. The SED was corrected for the interstellar extinction using E(B − V) = 0.15 mag (see Figure 5) and the average Galactic reddening law from Savage & Mathis (1979).

It was assumed that the IR excess radiation at wavelengths λ ≳ 1µm originates from a circumbinary dusty disk. The flux density in this case is given by

$$f_{\nu, \text{disk}} = d^{-2} \int_{R_{\text{in}}}^{R_{\text{out}}} B_{\nu}(T_r) (1 - \exp(-\tau_{\nu,r})) 2\pi r dr, \quad (1)$$

where R<sub>in</sub> and R<sub>out</sub> are the disk inner and outer radii, respectively, d is the distance to the star, and B<sub>ν</sub>(T) is the Plank function. The computations were done for a distance of 650 pc. Optical depth of the disk material τ<sub>ν,r</sub> is the product of a wavelength-dependent disk opacity, κ<sub>ν</sub>, and the radial surface density distribution, Σ<sub>r</sub>. To calculate the disk opacity, we used the Mie theory and considered spherical grains composed of astronomical silicates with a density of 2.5 g cm⁻³, sizes of a = 0.01–100 µm (see Boehler et al. 2013, for a detailed description of the grain emissivity determination), and a size distribution ∝ a⁻³.⁵ typical of interstellar grains (Mathis et al. 1977). We assumed that the disk has a gas to dust mass ratio of 100 (Bohlin et al. 1978). Contributions from both gas and dust are taken into account in calculating the opacity of the circumstellar material.

The temperature T<sub>r</sub>, vertical height H<sub>r</sub> and surface density Σ<sub>r</sub> distributions were parameterized as power laws of the disk inner radius:

$$T_r = T_{\text{in}} \left( \frac{r}{R_{\text{in}}} \right)^{-\bar{q}}, \quad (2)$$

$$H_r = H_{\text{in}} \left( \frac{r}{R_{\text{in}}} \right)^{\frac{3-q}{2}}, \quad (3)$$

$$\Sigma_r = \Sigma_{\text{in}} \left( \frac{r}{R_{\text{in}}} \right)^{-p}, \quad (4)$$

where T<sub>in</sub> is the disk temperature at R<sub>in</sub> that we calculate using the radiative equilibrium equation. The index q̄ and the disk vertical height at the inner edge H<sub>in</sub> are taken to be free parameters. The index p describes the surface density variation with distance from the star. Σ<sub>in</sub> is a surface density at R<sub>in</sub>. The total mass of the disk is related to the disk size and surface density as follows:

$$M_d = \int_{R_{\text{in}}}^{R_{\text{out}}} 2\pi r \Sigma_r dr. \quad (5)$$

Conversely, one can express Σ<sub>in</sub> as a function of the disk mass:

$$\Sigma_{\text{in}} = \frac{M_d R_{\text{in}}^{-p} (2-p)}{2\pi (R_{\text{out}}^{2-p} - R_{\text{in}}^{2-p})}. \quad (6)$$

Geometry of a disk tilted with respect to the line of sight is accounted for by the approach described in Zakhozhay et al. (2015). It assumes that the disk edges at the inner and outer
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Variations of the Hα line profile parameters folded with the orbital phase. Upper left panel: peak intensity ratio. Lower left panel: red emission peak. Intensity scale is given in continuum units. Dotted lines connect the measurements taken during the same orbital cycle: 2014—circles, 2015—squares, 2017—upward triangles, and 2018—downward triangles. Only the cycles with a representative phase coverage are shown. Upper right panel: the Hα line emission-peak separation. Lower right panel: the Hα line red emission peak position. RVs are shown in kilometers per second.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Table 3**
Parameters of the Dusty Disk</figcaption></figure>

<table><thead><tr><td>Parameter</td><td>Best-fit</td><td>Range</td><td>Increment</td></tr></thead><tbody><tr><td>R<sub>out</sub></td><td>60 au</td><td>20 &divide; 850 au</td><td>10 au</td></tr><tr><td>q&#773;</td><td>0.75</td><td>0.00 &divide; 0.85</td><td>0.05</td></tr><tr><td>p</td><td>1.5</td><td>0.0 &divide; 2.0</td><td>0.3</td></tr><tr><td>M<sub>d</sub></td><td>5 10<sup>-5</sup> M<sub>2</sub></td><td>10<sup>-5</sup> &divide; 10<sup>-1</sup>M<sub>2</sub></td><td>10%</td></tr><tr><td>H<sub>in</sub></td><td>0.08 R<sub>in</sub></td><td>10<sup>-3</sup> &divide; 10<sup>-1</sup>R<sub>in</sub></td><td>10%</td></tr></tbody></table>

**Note.** The mass of the mass gainer, $M_2$, is listed in Table 4. The uncertainties of the best-fit parameters are discussed in the text.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Part of the TCO spectrum of 3 Pup showing the diffuse interstellar bands (DIBs) that were used to estimate the amount of interstellar reddening. Intensities and wavelengths are in the same units as in Figure 1.</figcaption></figure>

radii have a cutoff flat geometry and emit as blackbodies with constant temperatures derived from the radiative equilibrium equation and Equation (2) at the inner and outer edges, respectively.

Using the modeling approach described above and tested on the SED of the young star IRAS 22150 + 6109 (Zakhozhay et al. 2018), a grid of SEDs for the disk with different parameters was calculated. The ranges and increments of the modeling parameters are shown in Table 3.

The best fit for the SED was found by minimizing:

$$\chi^2 = \sum_{i=1}^{n} \left( \frac{F_{\text{obs},i} - F_{\text{mod},i}}{F_{\text{obs},i}} \right)^2, \quad (7)$$

where $F_{\text{obs},i}$ and $F_{\text{mod},i}$ are the observed and modeled fluxes (at the corresponding wavelength) respectively. We normalize all the differences to the observed fluxes to account for their large

range (~$10^6$, see Figure 6). We compared the modeled and observed SED only at wavelengths >1 µm, where the disk dominates the object's radiation. At all tilt angles used in the fitting process, there was no attenuation of the A-supergiant by the disk. This is justified by a very good agreement of the dereddened observed UV fluxes and the intrinsic model SED for the chosen $T_{\text{eff}}$ of the star as well as by a much poorer agreement of the observed and model fluxes in the near-IR region at a very high tilt angle of 80°.

We fixed $R_{\text{in}}$ at 4 au, because this is consistent with the results of the interferometric observations (Meilland et al. 2010; Millour et al. 2011). This value is in agreement with the minimum possible radius, which is the dust sublimation radius (3.88 au), assuming that the dust sublimation temperature is 1500 K. The best-fit ($\chi^2$= 1.6, Figure 6) model parameters are listed in the second column of Table 3. Three of them were constrained reasonably well: $H_{\text{in}} = 0.08_{-0.05}^{+0.02} R_{\text{in}}$, $\tilde{q} = 0.75_{-0.15}^{+0.05}$, and $p = 1.5_{-0.5}^{+0.3}$. The error bars correspond to $1\sigma$ uncertainties calculated from the $\Delta\chi^2$ confidence statistics.

Since the disk is optically thin, which is supported by the presence of the silicate emission features at 10 and 18 µm (see the upper inset in Figure 6), the theoretical SED is not very sensitive to the tilt angle variation in a range of 40°–60° and still results in a reasonably good fit to the observed one. The fits
<figure><img src="image_9.png" /></figure>

**Figure 6.** Best-fit disk model (dashed line) is shown in comparison with the interstellar reddening corrected observed SED of 3 Pup composed from the mentioned photometric data (filled circles) and IRAS low resolution spectrum (open circles in the upper inset). The lower inset shows the agreement of the observed and model UV fluxes. The uncertainties of the observational data points do not exceed the symbol size except for those at $\lambda = 11.6, 12.8$, and $450~\mu$m, whose errors are ~30%. The photospheric radiation of the binary is shown with the solid line. It is represented by the scaled sum of the emergent fluxes from the mass gainer with $T_{\text{eff}} = 8500$ K and log g = 2.0 and the mass donor with $T_{\text{eff}} = 50,000$ K and log g = 5.0 from a model grid by Castelli & Kurucz (2003) for the solar element abundances. The component contributions are set by parameters from Table 4. The flux from the mass donor becomes dominant at wavelengths shorter than ~1000 Å. The dotted line shows a model SED of the mass gainer alone.

get worse outside of this angle range. Therefore, our results on the tilt angle are in good agreement with those from the previous modeling: $60^{\circ} \pm 10^{\circ}$ by Meilland et al. (2010) and $38^{\circ}$ Millour et al. (2011).

Similarly, the model does not constrain the disk outer radius, $R_{\text{out}}$, which determines the flux at the longest wavelengths. The minimum value of $R_{\text{out}}$ that gives a reasonably good agreement with the observed submillimeter flux is 50 au. Increasing it to 120 au only weakly contributes to the long-wavelength flux given the temperature decrease with distance from the star and a low optical depth of the dust. This uncertainty affects the total disk mass, which becomes an order of magnitude larger than the best-fit value presented in Table 3 at $R_{\text{out}} \sim 120$ au. Nevertheless, even with these problems the total mass of the circumstellar material (dust and gas) remains below $0.01~M_{\odot}$ and represents a very small fraction of the total mass of the entire system.

## 4. Discussion

Our orbital solution clearly suggests a stable and unique orbital period of 137.4 days. This is the same as suggested in the very first study of the RV variations by Johnson & Neubauer (1946) and does not support later results by Plets et al. (1995). The RV semiamplitude is small and, along with the orbital period, results in a small mass function (see Section 3.1). The mass function depends on three parameters: masses of both components and the tilt angle of the system rotational axis with respect to the line of sight:

$$f(m) = \frac{M_1 \times \sin^3 i}{(1 + q^2)}, \quad \text{where } M_1 \text{ is the mass of the currently lower-mass component, } q = M_2/M_1 \text{ is the component's mass ratio, } M_2 \text{ is the mass of the A-supergiant, and } i \text{ is the tilt angle.}$$

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

**Figure 7.** Hertzsprung-Russell diagram with positions of FS CMa objects with known fundamental parameters (see Miroshnichenko 2007). Evolutionary tracks for single rotating stars are taken from Ekström et al. (2012) and shown by solid lines along with the main-sequence location. Numbers by the tracks indicate initial masses in solar units. Recently identified binary systems are marked with their IDs.

In order to estimate the mass of the A-supergiant, its fundamental parameters need to be determined. The spectral type has varied from B8 I to A3 II in different papers summarized in SIMBAD, while the effective temperature was reported to be from 8250 K (Aret et al. 2016) to 10,000 K (Stee et al. 2004). Comparison with optical spectra of a number of A-type supergiants done by Chentsov et al. (2010) and in this study (see Figure 1) favors $T_{\text{eff}} = 8500$ K with an uncertainty on the order of 500 K. It is also consistent with the UV and visual SED corrected for the interstellar reddening $E(B-V) = 0.15 \pm 0.02$ mag, which was derived from spectral features, such as the 2175 Å band (Rovero & Ringuelet 1994) and diffuse interstellar bands (see Figure 5), whose EW were converted into $E(B-V)$ using a calibration by Herbig (1993).

The interstellar extinction calculated from the derived reddening, $A_V = 3.1 \times 0.15 = 0.47$ mag, the absolute visual magnitude of $M_V = -5.5 \pm 0.3$ mag from Chentsov et al. (2010), and the average visual magnitude ($V = 3.96 \pm 0.03$ mag) imply a distance of $630 \pm 85$ pc. The latter within the measurement uncertainty coincides with the one ($633_{-107}^{+160}$ pc) calculated from the GAIA parallax ($1.58 \pm 0.32$ mas, Gaia Collaboration et al. 2018). A very similar distance (650 pc) was used in all the papers published after 2004 quoted here.

The mass of the A-supergiant can be estimated from the evolutionary tracks for single stars, as no signs of the current mass transfer are seen in both photometric (no large brightness variations) and spectroscopic data (weak emission-line spectrum with small and gradual line profile variations). The adopted values of $M_V$ and $T_{\text{eff}}$ imply a luminosity of log $L/L_{\odot} = 4.1 \pm 0.1$. Comparison of these parameters with the evolutionary tracks for rotating single stars with initial atmospheric solar abundances (Ekström et al. 2012) suggests two possible values of the A-supergiant mass: $9.5 \pm 0.5~M_{\odot}$ for the pre-red-supergiant phase and $8.8 \pm 0.5~M_{\odot}$ for the post-red-supergiant phase (see Figure 7).

These mass estimates along with the orbital solution can be used to derive the mass of the secondary component. The dusty disk tilt angle derived from our SED modeling is consistent with the interferometric results (Meilland et al. 2010; Millour et al. 2011). Assuming that the binary orbit is tilted with respect
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Mass function of the 3 Pup binary system. $M_1$ and $M_2$ are masses of the donor and gainer in solar units, respectively. The darker-shaded area represents the relationship between the masses for the orbital tilt angles of 40° and 60° and the best-fit mass function ($f(m) = 1.81 \times 10^{-3} M_{\odot}$). The lighter-shaded area shows how the relationship changes for $f(m) \pm 1\sigma$ within the same interval of tilt angles. Vertical dashed lines show suggested limits for the mass of the donor based on the probable range of tilt angles. Horizontal dashed lines show the limits on the mass of the gainer from Table 4.</figcaption></figure>

to the line of sight at the same angle as the circumbinary disk ($i = 40°–60°$ from both interferometric and IR-excess modeling data), the secondary component should have a mass of $M_1 \sim 0.8 M_{\odot}$ (see Figure 8).

There can be two evolutionary scenarios that lead to the current masses of the stars in this system. One scenario assumes the components' evolution as single stars with no mass-exchange, while the other one involves mass-transfer between the components.

In the former scenario, the lower-mass component would be a solar-like star, which is much fainter than the A-supergiant, still located at the main sequence and not contributing to the observed spectrum. However, it seems to be less likely, because the A-supergiant is not massive enough to develop a strong stellar wind, which would supply a sufficient amount of material for the observed gaseous and dusty disk. Spectroscopic observations of many supergiants typically show either PCyg type profiles of Balmer lines (in the most luminous ones with the luminosity type Ia) or no line emission at all (in less luminous ones). Even the most luminous ones do not show the presence of dust in their circumstellar environments (e.g., Verdugo et al. 1999).

This reasoning applies to the earlier possible evolutionary phase of the A-supergiant (pre-red-supergiant). The possibility that it is currently a post-red-supergiant is very unlikely. When it reaches the lower-temperature end of the evolutionary track, its radius becomes so large that the system will go through a common envelope phase given the orbital period. After this phase the element abundances of the brighter component will be significantly altered (e.g., the hydrogen abundance becomes much lower, see Vanbeveren et al. 1998), while the observed absorption spectrum of 3 Pup is typical for a pre-red-supergiant phase star (see Figure 1). Therefore, we favor the scenario of the system evolution with mass-transfer between the components.

The derived components' masses and the presence of the circumstellar material are qualitatively consistent with the evolutionary models of binary systems that undergo mass

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Theoretical evolutionary tracks of a $6.0 + 3.6 M_{\odot}$ binary system with an initial orbital period of 5 days that we calculated using the code described in van Rensbergen et al. (2008). The dashed line shows the track for the mass donor, while the solid line shows that of the mass gainer. The numbers near the tracks indicate initial and final masses of the components in solar units, while the letters indicate the beginning and ending of the Roche overflow phase. The black circle shows the current position of the mass gainer of the 3 Pup binary system. The inset shows the mass transfer process in relative time since the beginning of the Roche lobe overflow phase.</figcaption></figure>

transfer calculated by van Rensbergen et al. (2008). A strong mass transfer during the Roche lobe overfilling time by a more massive and more evolved donor resulted in transferring most of its mass to the gainer as well as losing a fraction of the transferred mass to the circumstellar medium and beyond. Some models of van Rensbergen et al. (2008) show that the orbital period can reach several months after the end of the mass transfer phase, and that even the mass donor may get down to such a low mass as $0.6–1.0 M_{\odot}$.

Since there was no close match to the system parameters in the published model grid, we extended it using the same calculation approach and found one that has the final period (134 days) and luminosity of the more massive component close to those derived above. The adopted evolutionary track has initial masses of $6.0 + 3.6 M_{\odot}$ and an initial orbital period of 5.0 days using the same approach. It assumes a conservative mass transfer, since the mass of the circumstellar material that we found in our SED modeling (Section 3.3) is very small (see Figure 9). The Roche lobe overflow phase (partially shown in the inset of Figure 9) lasts for 0.4 Myr. The system is currently 9 Myr after this phase, $\sim 10\%$ of its total evolutionary time. These results along with a lack of studies of more distant B- and A-type supergiants may explain a rare detection of such objects.

The evolutionary track for the mass gainer gets close to the current position of the A-supergiant with a mass of $8.8 M_{\odot}$. After the end of the mass transfer phase, which lasts for $\sim 10^5$ yr, the mass gainer evolves as a single star with nearly the final mass (with a small mass loss through the stellar wind). The mass donor has a final mass of $0.8 M_{\odot}$, a $T_{\text{eff}} \sim 50,000$ K, and a bolometric luminosity an order of magnitude smaller than that of the mass gainer. With these parameters, the hot He-rich dwarf would noticeably contribute to the binary's SED only at wavelengths below $1000 \AA$ (see Figure 6). The spectrum of 3 Pup has not been observed in this far-UV region yet.

The difference in the single and binary evolutionary mass of the gainer is not large ($\sim 10\%$) and can be attributed to
<table><caption>**Table 4**
Parameters of the System Components</caption><thead><tr><th>Parameter</th><th>Gainer</th><th>Donor</th></tr></thead><tbody><tr><td>T<sub>eff</sub></td><td>8500 ± 500 K</td><td>~50,000 K</td></tr><tr><td>R/R<sub>☉</sub></td><td>54 ± 7</td><td>~0.3</td></tr><tr><td>log g</td><td>1.9 ± 0.1</td><td>~5.0</td></tr><tr><td>M/M<sub>☉</sub></td><td>8.8 ± 0.5</td><td>0.75 ± 0.25</td></tr></tbody></table>

Note. The parameters are based on our analysis of the spectrum of the A-supergiant and the evolutionary model of the binary (Figure 9). T<sub>eff</sub>, R, log g, and M are the current effective temperature, radius, surface gravity, and mass of each of the system components.

differences in the evolutionary codes of Ekström et al. (2012) and van Rensbergen et al. (2008) as well as to a different internal structure of the star that is set by the evolutionary path. Since the object’s binarity was revealed and confirmed, we adopt the binary evolutionary mass of $M_2 = 8.8 \pm 0.5 M_☉$ for the gainer. The mass uncertainty is based on that of the luminosity and should be similar to that estimated from the single star evolutionary tracks. The gainer’s mass leads to a slightly altered estimate for the donor’s mass ($M_1 = 0.75 \pm 0.25 M_☉$, see Figure 8). The final fundamental parameters of the system components are shown in Table 4.

Spectroscopically, the mass gainer in the 3 Pup binary system is a supergiant and has always been considered to belong to the B[e] supergiants subgroup. All other objects of the latter subgroup are typically much more massive with an average luminosity of an order of magnitude larger than that of the most luminous FS CMa objects (e.g., Miroshnichenko 2007). The position of this component in the Hertzsprung-Russell diagram (see Figure 7) implies that it is currently more consistent with the system classification as an FS CMa object. This subgroup of objects with the B[e] phenomenon was defined by Miroshnichenko (2007), and their properties were interpreted as due to consequences of the abovementioned nonconservative binary evolution. Most properties of 3 Pup (except for the strength of the emission-line spectrum) match those of the FS CMa group. Thus, we suggest a reclassification of the 3 Pup binary system into an FS CMa object.

The orbital solution and the components’ masses result in an orbital semimajor axis of $a = 1.11 \pm 0.03$ au or $4.4^{+0.7}_{-0.5} R_2$, where $R_2$ is the radius of the A-supergiant. According to Eggleton (1983), the Roche lobe of the mass gainer has a radius of $0.65 a$ for the current components’ mass ratio. Therefore, the mass gainer is well confined within its Roche lobe, and its gaseous disk is limited by the Roche lobe size to $\sim 2.9 R_2$.

The derived mass and radius of the mass gainer imply a critical rotation velocity $v_{crit} = 436.8 \sqrt{(M_2/M_☉)/(R_2/R_☉)} = 178^{+18}_{-16}$ km s<sup>-1</sup>. Our measurements with the Fourier transform method show that the projected rotational velocity of the A-supergiant is $v_{rot} = 35 \pm 5$ km s<sup>-1</sup>. This result indicates a slower rotation compared to the previously reported data (e.g., $50 \pm 5$ km s<sup>-1</sup>; Plets et al. 1995). Thus, one should not expect an enhanced mass loss from it. Assuming that the Keplerian gaseous disk extends all the way to the star and its density does not significantly drop outward, the separation of the line emission peaks should be $v_{peak} = 2 v_{crit} \sin i \sim 273$ km s<sup>-1</sup> (see, Huang 1972). However, as seen in Figure 4, the average peak separation in the Hα emission line is $\sim 120$ km s<sup>-1</sup>. The latter value is expected at a distance of $\sim 6 R_2$, near the L2 point of the system. Along with a

stable position of the Hα line central depression, this result may be interpreted as evidence for a circumbinary location of the gaseous disk. This idea was suggested by Plets et al. (1995), but the location of the second star in the system was unknown at that time.

The typically weaker blueshifted peak in this line can then be explained by the presence of a stellar wind from the A-supergiant. The presence of the wind is expected for such a large and relatively massive star, but, as we mentioned above, a purely wind-driven Hα emission is usually observed only in more luminous A-type Ia supergiants (e.g., Verdugo et al. 1999). Nevertheless, even a weak wind might provide enough mass flux into the disk to distort the observed Hα line profile. Modeling of the spectral line profiles is a complicated task and will be attempted in a follow-up paper.

The Hα line profile features described in Section 3.2 can be due to the following mechanisms. The higher level of the central depression along with the presence of the weak bluest peak (see Figure 3) may be due to an increase in the mass-loss rate from the A-supergiant. This process eventually results in mixing the new material with the disk and, as a consequence, a lower contribution of the circumstellar gas to the RV region responsible for the central depression. The V/R dependence on the orbital phase may be due to a higher contribution of the stellar wind to the line profile, since a smaller fraction of the part of the wind directed toward the observer is attenuated by the mass donor near the superior conjunction orbital phase. The redshift of the red emission peak along with a greater peak separation (see Figure 4) do not contradict the above explanation.

## 5. Conclusions

A long-term monitoring of the brightest object exhibiting the B [e] phenomenon, 3 Pup, resulted in refining the orbit of this binary system and detecting variations of the Hα line profile, the strongest emission feature in the object’s optical spectrum. From the spectroscopic and GAIA parallaxes, which agree with one another, the derived orbital solution, and consideration of the binary system evolutionary models, we determined the current masses of the system components to be $M_2 = 8.8 \pm 0.5 M_☉$ and $M_1 = 0.75 \pm 0.25 M_☉$. With such a mass, the mass gainer better fits into the group of FS CMa objects rather than into that of B[e] supergiants. The mass donor may be a helium-rich subdwarf of a much lower luminosity, explaining why no contribution from such a high-temperature object is observed in the UV fluxes of 3 Pup.

An evolutionary model of a binary system with the components’ initial masses of 6.0 and 3.6 $M_☉$ and conservative mass transfer was calculated to explain the current fundamental parameters of the A-type component and the observed orbital period. Although it does not imply mass loss from the system, the small amount of circumstellar material estimated from the IR-excess modeling is not expected to noticeably change the components’ evolution.

We also suggest that the circumstellar gaseous disk is circumbinary. This follows from the separations of the double-peaked emission-line profiles. A stellar wind, which is expected from the A-supergiant, can be responsible for the lower strength of the blueshifted peak in the Hα line profile. The circumstellar dust is located further away from the stars than the circumstellar gas, whose outer parts are most likely neutral
due to a low temperature of the A-supergiant and a weak
contribution of the ionizing radiation from the mass donor.

Modeling the object's IR excess allowed us to put new
constraints on the amount and structure of its dusty disk. The
best disk model agrees with the interferometry data on the tilt
angle with respect to the line of sight. It also suggests that a
dust mass of ~$10^{-5}-10^{-4} M_{\odot}$ is currently present in the
circumbinary disk.

Continuation of spectroscopic monitoring is strongly
suggested to increase the phase coverage during one orbital
cycle to investigate the line profile variations in more detail and
try to detect the onset of the next strengthening of the central
depression in the Hα line profile. Measurements of the system
far-UV fluxes (at $\lambda \le 1000 \AA$) are needed to verify the nature
of the mass donor.

We are grateful to the anonymous referee for useful
suggestions that allowed us to improve the paper, especially
the data interpretation. A.M. and S.Zh. acknowledge support
from DGAPA/PAPIIT Projects IN 100617 and IN 102120. V.
G.K. thanks the Russian Foundation for Basic Research for
partial support (project 18-02-00029a). M.P. acknowledges
support from the Russian Foundation for Basic Research
project 18-02-00554. The work was partially carried out within
the framework of Project No. BR05236322 "Studies of
physical processes in extragalactic and galactic objects and
their subsystems" financed by the Ministry of Education and
Science of the Republic of Kazakhstan. This paper is partly
based on observations obtained with the Canada-France-
Hawaii Telescope (CFHT), which is operated by the National
Research Council of Canada, the Institut National des Sciences
de l'Univers of the Centre National de la Recherche
Scientifique de France, and the University of Hawaii as well
as on observations obtained at the 2.7 m Harlan J. Smith
telescope of the McDonald Observatory (Texas, USA), and the
2.1 m of the Observatorio Astronómico Nacional San Pedro
Martir (Baja California, México). The observations at the
Canada-France-Hawaii Telescope were performed with care
and respect from the summit of Maunakea, which is a
significant cultural and historic site. This research has made
use of the SIMBAD database, operated at CDS, Strasbourg,
France.

ORCID iDs

A. S. Miroshnichenko https://orcid.org/0000-0002-3833-1038

S. A. Khokhlov https://orcid.org/0000-0001-5163-508X

References

Adelman, S. J., & Albayrak, B. 1997, IBVS, 4541, 1

Aret, A., Kraus, M., & Šlechta, M. 2016, MNRAS, 456, 1424

Boehler, Y., Dutrey, A., Guilloteau, S., & Piétu, V. 2013, MNRAS, 431, 1573

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Castelli, F., & Kurucz, R. L. 2003, in IAU Symp. 210, Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray (San Francisco, CA: ASP), A20

Chentsov, E. L., Klochkova, V. G., & Miroshnichenko, A. S. 2010, AstBu, 65, 150

Donati, J.-F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658

Eggleton, P. P. 1983, ApJ, 268, 368

Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146

Gaia Collaboration, Brown, A. G. A., Vallenari, A., Prusti, T., et al. 2018, A&A, 616, A1

Herbig, G. H. 1993, ApJ, 407, 142

Huang, S.-S. 1972, ApJ, 171, 549

Johnson, H. L., & Neubauer, F. J. 1946, PASP, 58, 248

Jura, M., Webb, R. A., & Kahane, C. 2001, ApJL, 550, L71

Klochkova, V. G., Sendzikas, E. G., & Chentsov, E. L. 2015, AstBu, 70, 99

Kraus, M., Oksala, M. E., Cidale, L. S., et al. 2015, ApJL, 800, L20

Manset, N., & Donati, J.-F. 2003, Proc. SPIE, 4843, 425

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Meilland, A., Kanaan, S., Borges Fernandes, M., et al. 2010, A&A, 512, A73

Merrill, P. W. 1934, PASP, 46, 156

Millour, F., Meilland, A., Chesneau, O., et al. 2011, A&A, 526, A107

Miroshnichenko, A. S. 2007, ApJ, 667, 497

Monnier, J. D., Millan-Gabet, R., Tuthill, P. G., et al. 2004, ApJ, 605, 436

Panchuk, V. E., Klochkova, V. G., & Yushkin, M. V. 2017, ARep, 61, 820

Plets, H., Waelkens, C., & Trams, N. R. 1995, A&A, 293, 363

Rovero, A. C., & Ringuelet, A. E. 1994, MNRAS, 266, 203

Saio, H., Georgy, C., & Meynet, G. 2013, MNRAS, 433, 1246

Savage, B. D., & Mathis, J. S. 1979, ARA&A, 17, 73

Stee, P., Bittar, J., & Lopez, B. 2004, ApJ, 602, 978

Sterken, C., Manfroid, J., Anton, K., et al. 1993, A&AS, 102, 79

Swings, P. 1950, AnAp, 13, 114

Tull, R. G., MacQueen, P. J., Sneden, C., & Lambert, D. L. 1995, PASP, 107, 251

van Rensbergen, W., De Greve, J. P., De Loore, C., & Mennekens, N. 2008, A&A, 487, 1129

Vanbeveren, D., De Loore, C., & Van Rensbergen, W. 1998, A&ARv, 9, 63

Verdugo, E., Talavera, A., & Gómez de Castro, A. I. 1999, A&A, 346, 819

Zakhozhay, O. V., del Burgo, C., & Zakhozhay, V. A. 2015, AASP, 5, 33

Zakhozhay, O. V., Miroshnichenko, A. S., Kuratov, K. S., et al. 2018, MNRAS, 477, 977