The Transformative Journey of HD 93521

Most galactic massive stars are found close to their birthplace in the disk, but we do find some examples of O- and B-type stars far from the disk (Keenan 1992). These remote objects are probably runaway stars that were ejected from the disk by a supernova explosion in a binary system or by close gravitational encounters with binary stars in their natal star clusters (Blaauw 1961; Gies & Bolton 1986; Hoogerwerf et al. 2001). There are several known examples of runaway binary stars in the halo (e.g., IT Lib; Martin 2003; Wysocki et al. 2022) that may be the result of a supernova in a triple system (Gao et al. 2019) or the dynamical ejection of a close binary (Leonard & Duncan 1990). Comparisons of the computed time of flight from a starting point in the disk to the current position with an estimate of the evolutionary age generally show that the massive halo stars have had sufficiently long lifetimes to reach their locations (Silva & Napiwotzki 2011; Raddi et al. 2021). However, there are some mysterious cases that appear to be younger than their time of flight, in conflict with the idea that they formed in the disk (Perets 2009).

An important example of a case with a lifetime shorter than the time of flight is the high Galactic latitude star HD 93521 (O9.5 IIInn; Sota et al. 2011). This star is over 1 kpc from the Galactic plane, and it has served as a distant lantern to explore the interstellar medium along the line of sight (Gringel et al. 2000; Savage et al. 2001). The spectrum of HD 93521 displays very broad photospheric line profiles that indicate very fast rotation and therefore a close to equatorial-on inclination (projected rotational velocity ${v}_{e}\sin i=435$ km s⁻¹; Howarth & Smith 2001). Detailed spectral analyses show that the helium lines are unusually strong, indicating a He enrichment due to internal mixing or some other process (Howarth & Smith 2001; Rauw et al. 2012; Cazorla et al. 2017). The spectral features display systematic line-profile variations that are associated with at least two modes of nonradial pulsation (Fullerton et al. 1991; Howarth & Reid 1993; Howarth et al. 1998; Rauw et al. 2008; Rzaev & Panchuk 2008; Rauw & Nazé 2021). Howarth & Reid (1993) were the first to show that the evolutionary lifetime of HD 93521 is much shorter than the estimated time of flight from the disk, opening the possibility that it is a rare massive star formed outside the disk.

Here we examine the origin of HD 93521 through a calculation of the time of flight (Section 2) and a reassessment of its evolutionary age based upon the recent distance determination from Gaia EDR3 (Section 3). We confirm the marked discrepancy between the evolutionary and kinematical timescales, and we show how this difference may be explained by assuming that the object was ejected from the disk as a close binary star that subsequently merged to become a rapidly rotating, single star (Section 4).

We calculated the time of flight from an origin in the disk using a numerical integration of the star's motion in the Galaxy. We used a scheme described in earlier work (Boyajian et al. 2005) that is based upon a model for the Galactic potential from Dehnen & Binney (1998). The starting parameters for the calculation are the star's position, distance, proper motions, and radial velocity. The astrometric data were taken from the Gaia EDR3 catalog (Gaia Collaboration et al. 2016, 2021), and we used a distance from Gaia EDR3 of $d={1.246}_{-0.102}^{+0.136}$ kpc (Bailer-Jones et al. 2021). This distance represents a significant downward revision from the spectroscopically derived value of 2.2 kpc (Howarth & Smith 2001), and we explore the consequences of this change in the next section.

Table 1 lists estimates of the star's radial velocity from a number of investigators. The large range in results is due to the difficulty in measuring the positions of the broad and shallow absorption lines in the spectra that appear distorted through the influence of the nonradial pulsations (Rauw et al. 2012). We added to these results radial velocity measurements we made on two sets of spectra. The first is a set of five high S/N and high spectral resolving power spectra obtained with the ESPaDOnS instrument on the Canada–France–Hawaii Telescope (Donati et al. 2006). The ESPaDOnS spectra were obtained from the PolarBase archive ³ (Petit et al. 2014), and were rebinned onto a heliocentric, $\mathrm{log}\lambda$ wavelength grid. We created a model spectrum template from the TLUSTY/SYNSPEC grid (Lanz & Hubeny 2003) on the same wavelength grid assuming stellar parameters given in the next section (Table 2 below). We then used the range lacking strong interstellar lines (3988–5886 Å) to cross-correlate each spectrum with the model, and the peak of the resulting cross-correlation function gave the radial velocity measurement. The average of these measurements is listed in Table 1. Finally, we applied the same kind of cross-correlation measurements to a set of Short Wavelength Prime camera, high dispersion, ultraviolet spectra from the archive of the International Ultraviolet Explorer (IUE) satellite. We selected those regions of the FUV spectrum that were free from strong wind features and interstellar lines to perform the cross correlation with a TLUSTY model spectrum, and the average of the resulting 139 radial velocity measurements appears in Table 1.

The columns of Table 1 list the source of the radial velocity measurements, the number of individual measurements, the mean date of the observations, the average radial velocity 〈V_r 〉, the external standard deviation between the measurements σ (E), and the average of the internal estimated errors for the individual measurements σ (I) (where available). A comparison of the latter two statistics shows that the external error is comparable to or slightly larger than the internal error as expected for small temporal variations that result from the influence of nonradial pulsations. Consequently, we assume the star is radial velocity constant within the measurement errors, and we formed an error-weighted mean from all the entries of Table 1, 〈V_r 〉 = −8.9 ± 2.9 km s⁻¹. The standard error of the mean in this estimate is based upon the number of degrees of freedom and the reduced ${\chi }_{\nu }^{2}$ that was renormalized to a minimum of unity. We used this radial velocity estimate for the galactic motion calculation.

The integration of the stellar motion backwards in time is illustrated in Figure 1 in Galactic (R, z) coordinates where R is the distance from Galactic Center and z is the distance from the plane. The star appears to have originated beyond the solar circle in a region occupied by the Sagittarius arm in quadrant four of the Galaxy. The star is now just beyond its maximum arc position above the plane and is beginning its return journey. The derived time of flight from the plane is 39 ± 3 Myr where the uncertainty is based upon the astrometric and radial velocity uncertainties. The ejection velocity from the disk was approximately 61 ± 3 km s⁻¹ relative to its local standard of rest following the galactic rotation curve. We also checked on systematic errors in the model results in two ways. By setting the starting position to z = ± 40 pc (the scale height for local OB associations; Bobylev & Bajkova 2016), the time of flight changed by ±1.2 Myr. Next we calculated the trajectory using another code galpy that uses a different description of the galactic potential (MWPotential2014; Bovy 2015). This integration of the motion arrived at the same estimate of the time of flight of 39 ± 3 Myr. Thus, the systematic uncertainties are comparable to the observational errors.

Zoom In Zoom Out Reset image size

Next we will compare the time of flight to the star's evolutionary age through inspection of the star's $({T}_{\mathrm{eff}},\mathrm{log}L/{L}_{\odot })$ position in the Hertzsprung–Russell (H–R) diagram relative to evolutionary track models for massive rotating stars (Brott et al. 2011). It is important to review the new constraints on the star's radius and luminosity that are implied by the recent, closer estimate of the star's distance. We begin by determining the star's angular diameter through a comparison of the star's spectral energy distribution with a model spectrum for a nonrotating star, and then we perform a more careful estimate of the model V-band monochromatic flux based on a rotating star model in order to estimate the stellar radius.

The spectral energy distribution of HD 93521 in the ultraviolet and optical part of the spectrum is well established thanks to its inclusion as a spectral flux standard in the Hubble Space Telescope (HST) CALSPEC ⁴ program (Bohlin et al. 2014, 2020). The source spectra are from IUE (1148–1680 Å) and HST/STIS (1680–10120 Å), and we rebinned these onto a $\mathrm{log}\lambda$ grid for a resolving power of R = 500. We then created a model SED with the same resolving power using fluxes from TLUSTY for assumed, hemisphere-averaged parameters of effective temperature T_eff = 30.46 kK and gravity $\mathrm{log}g=3.66$ (Rauw et al. 2012; Cazorla et al. 2017). Figure 2 shows the observed and model SEDs that agree well for fitting parameters of an angular diameter of θ = 0.0482 ± 0.0012 milliarcsec and a reddening of E(B − V) = 0.053 ± 0.007 mag (R_V = 3.1 assumed). Then, using the Gaia EDR3 distance, we find a mean radius of 6.5R_⊙ for HD 93521.

Zoom In Zoom Out Reset image size

We know, however, that the star is rapidly rotating, so its equatorial radius will be larger than its polar radius. Thus, in order to estimate the star's luminosity more accurately, we need to compare its observed flux with model predictions based on a realistic rotational model. We created such a model using a code that was originally developed to study the rapid rotator Regulus (McAlister et al. 2005) and that is described in detail by Shepard et al. (2022). This code simulates the photosphere of a rapidly rotating star by calculating the flux from each of the surface elements over the visible hemisphere. The rotational distortion is described by the Roche model, and the gravity darkening is assumed to follow the ω-model prescription from the work of Espinosa Lara & Rieutord (2011, 2013). The monochromatic specific intensities associated with the local temperature, gravity, and viewing angle of each surface element are derived from the TLUSTY/SYNSPEC models (Lanz & Hubeny 2003, 2007). The final hemisphere-integrated flux is then rescaled for the assumed distance and interstellar extinction (from E(B − V) = 0.053 mag derived above) to create a predicted flux that can be compared to the observed value.

Howarth & Smith (2001) applied a similar kind of model to fit the spectral lines of HD 93521 and derive the stellar properties. Because their model provided a successful fit of the observed spectrum, we decided to adopt many of their fitting parameters in our rotational model. However, a number of the model parameters needed revision. Howarth & Smith (2001) relied on atmosphere models that ignored the effects of line blanketing, and subsequent models that include line blanketing (like TLUSTY) usually arrive at a lower effective temperature in fits of stellar spectra. Since our model is based on specific intensity spectra from TLUSTY, it was important to revise downward the polar effective temperature T_p so that the hemisphere-average temperature matched estimates from recent, line-blanketed model fits of the spectrum of HD 93521 (Rauw et al. 2012; Cazorla et al. 2017). We also decided to use the physically motivated ω-model instead of the simpler von Zeipel gravity darkening law that was adopted by Howarth & Smith (2001). Finally, we needed to lower the estimated polar radius R_p in order to reduce the integrated flux for the closer distance from Gaia EDR3. This meant that we also needed to rescale the assumed mass M ∝ R_p in order to maintain the fraction of angular rotation to the critical value, Ω/Ω_c , derived by Howarth & Smith (2001). Note that this choice results in poorer agreement with their gravity estimates because of the different scaling, $M\propto {R}_{p}^{2}$, for constant gravity.

Our model fit was set to match the observed monochromatic flux at a wavelength of 5450 Å, F_λ = 5.85 × 10⁻¹² erg cm⁻² s⁻¹ Å⁻¹. The fitting parameters are listed in Table 2 and an image of the model star appears in Figure 3. Table 2 lists columns for the physical parameter and the values estimated from Howarth & Smith (2001) and this work. Howarth & Smith (2001) give representative uncertainties for the parameters that broadly apply to all three stars in their study, and we include them here for completeness. The first row indicates the gravity darkening law adopted in the analysis, and the next five rows give the basic geometrical properties of the model. These include the projected rotational velocity ${v}_{e}\sin i$, the inclination of the spin axis to the line of sight i, the equatorial velocity v_e , the critical velocity at which gravity and centrifugal acceleration balance v_c , and the ratio of the angular velocity to that for critical rotation Ω/Ω_c . We have adopted all of these five parameters directly from the work of Howarth & Smith (2001). The next row gives the polar radius R_p that we derive from the fit of the V-band flux and the adopted values of distance and extinction. The uncertainty for R_p is based upon the fractional errors in distance, angular size, and average temperature. The next two rows report the equatorial radius R_e (which is derived from R_p and the rotational geometry) and the mass M. The mass is rescaled from the Howarth & Smith (2001) estimate in proportion to R_p , so the fractional error remains the same in dex. The radius and mass are lower in our model as a result of the smaller radius implied by the closer distance from Gaia EDR3. The next three rows list the logarithm of gravity at the pole $\mathrm{log}{g}_{p}$, equator $\mathrm{log}{g}_{e}$, and the flux-weighted average over the visible hemisphere $\mathrm{log}g$(avg). The larger errors given for our gravity results reflect the large uncertainty in mass. The next four rows give the effective temperature for the pole T_p , equator T_e , the flux-weighted average over the visible hemisphere 〈T〉(avg),and the area-integrated average over the entire star 〈T〉(all). These values are all scaled from 〈T〉(avg) that is set to the observed value from recent spectroscopic analyses. The next row presents the luminosity integrated over the area of the entire surface $\mathrm{log}(L/{L}_{\odot })$, and its uncertainty results from the errors in R_p and 〈T〉(all). The penultimate row gives the number ratio y of He to H in the atmosphere. Note that our choice of He abundance was governed by the closest value available in a preconstructed grid of specific intensities (Shepard et al. 2022), and it plays no role in the fitting process. The final row lists the adopted distance d.

Zoom In Zoom Out Reset image size

We can now use the revised stellar parameters to estimate the evolutionary age within the framework of stellar evolution models for rotating massive stars developed by Brott et al. (2011). Figure 4 shows the evolutionary tracks for initial masses of 15 and 20M_⊙ that bracket our revised mass of 17M_⊙. Both of these models are based on starting equatorial velocities that are close to that of HD 93521 and that roughly maintain the same value during the H-core burning, main-sequence phase of evolution. Our derived estimates of $({T}_{\mathrm{eff}},\mathrm{log}L/{L}_{\odot })$ for HD 93521 are indicated by a diamond symbol in Figure 4, and its position is consistent with expectations for a slightly evolved star of mass 17M_⊙. Comparison with the evolutionary tracks implies an evolutionary age of 5 ± 2 Myr for HD 93521. This is far less than the time of flight of 39 ± 3 Myr (Section 2), so we confirm that there is a large discrepancy between the star's apparent youth and the long travel time if it was ejected from the disk.

Zoom In Zoom Out Reset image size

Inspection of Figure 1 shows that HD 93521 was located about 1.1 kpc above the plane at a time 5 Myr in past. We briefly consider whether it may have formed at such a great distance from the molecular clouds of the disk where massive star formation generally occurs. The halo of the Galaxy does contain regions of relatively higher density gas, particularly in radio-detected high-velocity clouds. However, searches for young stars in these clouds indicate that star formation is very rare and restricted to high density environments (Ivezić & Christodoulou 1997; Simon & Blitz 2002; Stark et al. 2015). Figure 1 shows that HD 93521 had only a small velocity in the direction normal to the disk at a time 5 Myr ago, so it carries no imprint of an infalling natal cloud. The M-cloud complex appears in the sky in the vicinity of HD 93521, but it is actually situated much farther away than the star (Danly et al. 1993). Rauw et al. (2012) found that the immediate region around HD 93521 is devoid of any nebulosity or low mass, X-ray emitting stars that often are found in regions of massive star formation. Consequently, there is no evidence to support the idea that HD 93521 was formed in the halo near its current location. We discuss a compelling alternative explanation in the next section.

Table 2. Stellar Parameters

Parameter	Howarth & Smith (2001)	This Paper
Gravity darkening	von Zeipel	ω-model
${v}_{e}\sin i$ (km s−1)	435 ± 20	435 ± 20
i (degrees)	90 ± 10	90 ± 10
ve (km s−1)	435 ± 20	435 ± 20
vc (km s−1)	⋯	595 ± 25
Ω/Ωc	${0.9}_{-0.10}^{+0.05}$	${0.9}_{-0.10}^{+0.05}$
Rp (R⊙)	9.7 ± (0.3 dex)	6.1 ± 0.3
Re (R⊙)	⋯	7.4 ± 0.4
M (M⊙)	27 ± (0.4 dex)	17 ± (0.4 dex)
$\mathrm{log}{g}_{p}$ (dex cgs)	3.9 ± 0.1	4.1 ± 0.4
$\mathrm{log}{g}_{e}$ (dex cgs)	3.5 ± 0.1	3.7 ± 0.4
$\mathrm{log}g$(avg) (dex cgs)	⋯	3.8 ± 0.4
Tp (kK)	38.0 ± 1.5	34.6 ± 1.2
Te (kK)	29.9 ± 1.5	28.7 ± 1.2
〈T〉(avg) (kK)	33.5 ± 1.5	30.5 ± 1.2
〈T〉(all) (kK)	⋯	31.0 ± 1.2
$\mathrm{log}(L/{L}_{\odot })$	5.1 ± 0.5	4.6 ± 0.1
y = N(He)/N(H)	0.18 ± 0.03	0.20
d (kpc)	2.2	${1.246}_{-0.102}^{+0.136}$

Download table as:  ASCIITypeset image

The basic solution to the timescale problem is based on the idea that the current stellar mass was originally distributed in two lower-mass stars with longer main-sequence lifetimes (Perets 2009; Rauw et al. 2012). Thus, we suggest that HD 93521 is the product of a stellar merger of two stars that were ejected as a close binary system from a starting location in the disk. We posit that the system began as a close (perhaps contact) binary that was ejected through dynamical interactions in a multiple stellar system or through the gravitational release from a neighboring massive star that exploded as a supernova. This system continued to move into the halo as a binary system until the more massive component reached the end of core H burning and engulfed its companion during expansion. If this merger reset the H-burning clock of the combined star, then the merger happened relatively recently, about 5 Myr ago. This line of reasoning implies that the primary component of the initial binary must have had a mass low enough to last through the difference between the time of flight and merger event, about 34 Myr. According to the moderate rotation rate tracks from Brott et al. (2011), a main-sequence lifetime of 34 Myr occurs for a star with an initial mass of about 8M_⊙. Thus, the binary may have started life as an 8M_⊙ plus 8M_⊙ pair that merged with little mass loss to arrive at a current mass of about 16M_⊙, which is slightly below our mass estimate of 17M_⊙.

The outcome of a merger event is a subject of active investigation (Glebbeek et al. 2013; Jiang et al. 2013; Schneider et al. 2019; Menon et al. 2021). Most studies predict that the merger product will be a fast rotator experiencing mass loss, and HD 93521 shows evidence of both stellar wind mass loss (Bjorkman et al. 1994) and Hα emission from mass loss into an equatorial disk (Rauw et al. 2008; Rzaev & Panchuk 2008). Binary mass transfer may transport nuclear processed gas into the atmosphere of the mass gainer star, producing an enrichment in helium. Spectroscopic investigations indicate that there is a significant He enrichment in the atmosphere of HD 93521 (Howarth & Smith 2001; Rauw et al. 2012). Schneider et al. (2019) recently presented a pertinent numerical simulation of the merger of a 9M_⊙ and 8M_⊙ pair that forms a rapidly rotating star of 17M_⊙. They use a 3D magneto-hydrodynamical simulation to follow the gas flows created during the merger and the magnetic fields that are generated in the process. These simulations suggest the merger product will maintain a magnetic field, which may eventually cause the star to spin down. This prediction may be consistent with the 200–300 G longitudinal magnetic field of HD 93521 detected by Hubrig et al. (2013).

Finally the properties of the nonradial pulsations may offer another clue about the binary origin of HD 93521. Lee & Saio (2020) argue that pulsations can be excited in the envelope of a rapid rotator if the core rotates at a slightly faster rate than the envelope, and the superperiod mP, where m is the azimuthal order and P is the pulsation period, should be similar to the core rotation period. The superperiods found by Howarth et al. (1998) are 17.7 and 17.4 hr for the P = 1.77 (m = 10) and P = 2.90 hr (m = 6) modes, respectively, and these are less than the estimated rotation period of the photosphere of 20.7 hr (based on our solution in Table 2). Thus, HD 93521 may represent a case where a slightly faster rotating core drives envelope pulsations. We can compare this core rotational period to the expected orbital period of the premerger binary of ≈ 18 hr for an assumed semimajor axis of a = R₁ + R₂ = 9R_⊙. This short orbital period is comparable to that found in massive close contact binaries (Yang et al. 2019). We expect that both stars would attain synchronous rotation in the premerger binary because of the strong tidal forces, so the spin periods of the precursor stars would be the same as the orbital period. Thus, if we identify the pulsational superperiod as the core rotation period, then it appears that the core is rotating with a period that is about the same as the rotation and orbital period of the premerger stars.

The observed properties of HD 93521 all appear to agree with expectations for a merger product. The star appears to be too young compared to its time of flight from the Galactic disk because it was rejuvenated through the stellar merger of the binary components. Rejuvenation by a binary interaction may provide a common explanation for other such massive stars in the halo (Perets 2009). For example, Wysocki et al. (2022) recently found that the halo eclipsing binary IT Librae is in fact a post-mass-transfer system with a rejuvenated and massive primary star. Investigations of such systems will provide important clues about the properties of post-mass-transfer and merger systems that are key to understanding their final supernova progeny (Chrimes et al. 2020).

We are grateful for helpful comments from an anonymous referee. This work is supported by the National Science Foundation under grant No. AST-1908026. Institutional support has been provided from the GSU College of Arts and Sciences and the GSU Office of the Vice President for Research and Economic Development. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Facilities: IUE - International Ultraviolet Explorer, HST (STIS) - , CFHT (ESPaDOnS). -

Software: galpy (Bovy 2015), TLUSTY/SYNSPEC (Hubeny & Lanz 2017).