Kinematics of the Central Stars Powering Bowshock Nebulae and the Large Multiplicity Fraction of Runaway OB Stars

OB stars powering stellar bowshock nebulae (SBNe) have been presumed to have large peculiar velocities. We measured peculiar velocities of SBN central stars to assess their kinematics relative to the general O star population using $Gaia$ EDR3 data for 267 SBN central stars and a sample of 455 Galactic O stars to derive projected velocities $v_{\rm 2D}$. For a subset of each sample we obtained new optical spectroscopy to measure radial velocities and identify multiple-star systems. We find a minimum multiplicity fraction of 36$\pm$6% among SBN central stars, consistent with $>$28% among runaway Galactic O stars. The large multiplicity fraction among runaways implicates very efficient dynamical ejection rather than binary-supernova origins. The median $v_{\rm 2D}$ of SBN central stars is $v_{\rm 2D}$=14.6 km s$^{-1}$, larger than the median $v_{\rm 2D}$=11.4 km s$^{-1}$ for non-bowshock O stars. Central stars of SBNe have a runaway ($v_{\rm 2D}$$>$25 km s$^{-1}$) fraction of 24$^{+9}_{-7}$%, consistent with the 22$^{+3}_{-3}$% for control-sample O stars. Most (76%) of SBNe central stars are not runaways. Our analysis of alignment ($\Delta_{\rm PA}$) between the nebular morphological and $v_{\rm 2D}$ kinematic position angles reveals two populations: a highly aligned ($\sigma_{PA}$=25$^\circ$) population that includes stars with the largest $v_{\rm 2D}$ (31% of the sample) and a random (non-aligned) population (69% of the sample). SBNe that lie within or near HII regions comprise a larger fraction of this latter component than SBNe in isolated environments, implicating localized ISM flows as a factor shaping their orientations and morphologies. We outline a new conceptual approach to computing the Solar LSR motion, yielding [U$_\odot$, V$_\odot$, W$_\odot$]= [5.5, 7.5, 4.5] km s$^{-1}$.


Bowshock Nebulae and their Central Stars
Stellar bow shock nebulae (SBNe) constitute a distinctive class of arcuate interstellar nebulae flanking a central star traveling supersonically through the interstellar medium, first identified in association with ζ Oph and LL Ori (Gull & Sofia 1979). These nebulae contain swept-up stellar wind material and trapped ISM at a circumstellar distance where the momentum flux of the ISM balances that of the stellar wind (Wilkin 1996(Wilkin , 2000. The compressed and heated material of the nebula is most easily seen in the infrared continuum, though some SBNe have detectable emission lines including Hα or O III (Gull & Sofia 1979;Brown & Bomans 2005). As more SBNe have been discovered with advances in angular resolution and sky coverage of infrared surveys, these nebulae have generated increasing interest as laboratories to explore the properties of both the central star and the ISM interactions forming the nebular structure. SBNe and SBN-like objects surround a wide variety of stars including red supergiants (Noriega-Crespo et al. 1997;Gvaramadze et al. 2014), pulsars (Wang et al. 2013), externally ionized proplyds (Bally et al. 1998), and even an A dwarf (δ Vel, Gáspár et al. 2008). Chick et al. (2020, hereafter Paper I) conducted a spectroscopic survey of 84 cataloged SBN-candidates (Kobulnicky et al. 2016, hereafter K16) and determined that >95% were O and early B stars, in agreement with theoretical expectations (Weaver et al. 1977;Mac Low et al. 1991;Comeron & Kaper 1998).
Hydrodynamical simulations indicate that the density of the ambient ISM and the velocity differential between the ISM and the central star are both significant factors in the creation of an observable SBN (Comeron & Kaper 1998;Meyer et al. 2014Meyer et al. , 2016. Observable SBNe likely require ambient interstellar densities of 1-5 cm −3 (Kobulnicky et al. 2018;Henney & Arthur 2019), making it less likely that SBNe can form or be detected inside H II regions that efficiently heat and expel the ISM material (Hills 1980;Goodwin & Bastian 2006;Dale et al. 2012) or in other rarefied environments. The star-ISM velocity differential may be provided by the motion of the star (van Buren & McCray 1988;van Buren et al. 1990;Mac Low et al. 1991;van Buren & Mac Low 1992) or by a bulk ISM flow impinging upon a stationary star (i.e., "in-situ" bowshocks; Povich et al. 2008), and perhaps a combination of both. Alternatively, the stellar wind momentum flux may be supplemented (or even dominated) by stellar radiation pressure (van Buren & McCray 1988), particularly around late-O stars in the "weak wind regime" (Ochsendorf et al. 2014a,b;Henney & Arthur 2019). The diverse range of interaction types has given range to an equally diverse nomenclature of SBN subtypes, one still in the process of being standardized. Some authors (Henney & Arthur 2019) have suggested three major subtypes: traditional "bow shocks" where a thick shell of gas and dust forms about a dust-free shocked wind-supported (or radiation-supported) inner layer, "bow waves" where a thin shell of unshocked gas and dust forms due to the star's radiation pressure acting on the dust and a gradual deceleration by gas-dust drag, and "dust waves" where the gas and dust decouple and a thin dusty shell forms supported by radiation pressure but ISM gas continues to flow. Variations and hybrids of each of these models have been proposed as well (van Marle et al. 2011). Throughout this work, we will refer to all of these arc-shaped nebulae as SBNe, acknowledging that this sample may contain some or all of these physically distinct subtypes. Henney & Arthur (2019) analyzed the 20 SBNe from Kobulnicky et al. (2018) and found 15 to be bow shocks, two likely bow waves, and three ambiguous objects.
The presumption that SBNe require supersonic (but not hypersonic) star-ISM velocity differentials and the need to form in denser regions of the ISM than their natal clusters suggests that "runaway" stars (Blaauw 1961) are good candidates to produce SBNe. Runaways exhibit large peculiar velocities, having been ejected from their natal clusters via N-body dynamical interactions (Poveda et al. 1967;Spitzer & Mathieu 1980) or by supernovae in binary systems (Zwicky 1957;Blaauw 1961). Indeed, K16 found 70% of their SBN central stars were in "isolated" environments away from H II regions. A minority of sources (7%) were co-located on the sky with H II regions. Although the threshold for runaway velocities varies by author, most studies agree on velocities between 28 km s −1 (Tetzlaff et al. 2011) and 40 km s −1 (Blaauw 1961), with most modern works adopting a three-dimensional peculiar velocity threshold of 30 km s −1 (Gies & Bolton 1986). Additionally, most runaways are OBA-type stars (Tetzlaff et al. 2011;de la Fuente Marcos & de la Fuente Marcos 2018), though this may be a selection effect due to their large luminosities. While runaways provide a seemingly natural solution to both the ISM density and velocity differential requirements to form an SBN, it is clear that not all central stars supporting SBNe are runaways and not all runaways support observable SBNe. Some studies have found SBNe that are aligned with ISM flows emanating from H II regions, suggesting the presence of an in-situ SBN (Gull & Sofia 1979;Povich et al. 2008). Similarly, Bodensteiner et al. (2018) found candidate SBNe that show significant misalignment between their nebular morphological axis and their stellar proper motion vectors. Peri et al. (2012Peri et al. ( , 2015 searched archival infrared images of known OB runaway stars and found that only 15% produced SBNe observable in the infrared continuum. Based on these works, it appears that detectable SBNe require a specialized set of circumstances. Being a runaway OB star alone is not sufficient.

Runaway Stars and their Origins
Observations indicate that runaway and "walkaway" stars (peculiar velocities <30 km s −1 ) comprise the majority of field OB stars-those lying outside known H II regions and young clusters. de Wit et al. (2005) found that only 4% of field stars could not be associated known associations and star forming regions, implying that very few form in isolation. While some additional OB stars may have formed in isolated environments (Oey et al. 2018), field and cluster OB stars are statistically indistinguishable, except that runaway O stars have later spectral types than those in clusters (van den Bergh 2004). It is likely that isolated formation of O stars, if possible, is extremely uncommon.
Kinematic investigations O stars have yielded significantly different runaway fractions depending on the particular observational selection criteria (e.g., definition of runaway velocity threshold) or theoretical modeling approach. Even studies that draw samples of stars from the same catalogs can generate disparate results. For example, based upon the radial velocities of stars from the catalogs of Rubin et al. (1962) and Cruz-González et al. (1974), Stone (1979) used a 40 km s −1 velocity threshold to find an OB runaway fraction of 49±10%, whereas Gies & Bolton (1986) used a 30 km s −1 velocity threshold to find a runaway fraction of 10-25% for O stars and 2% for early B stars. The majority of studies identify runaways based upon peculiar radial velocities or peculiar transverse velocities alone. Few studies incorporate all three velocity components. Assuming the distribution of stellar velocities is isotropic, studies based upon the two components of proper motion instead alone will underestimate the true runaway fraction by 22% (i.e., 3/2; Lamb et al. 2016). By a similar reasoning, studies based upon radial velocity alone may underestimate the true fraction by as much as 73% (i.e., 3/1).
In the seminal work on runaway O stars, Blaauw (1961) adopted a 3D runaway threshold of 40 km s −1 to find a runaway fraction of 1.5% for early-B stars, 2.5% for B0 stars, and 21% for late-O stars. Stone (1979) reanalyzed archival Yerkes Observatory plates to measure the 2D peculiar velocities of 47 O4-B2 stars. Using a runaway threshold of 40 km s −1 they found an OB runaway fraction of 49±10%, though Gies & Bolton (1986) note that this sample may have been biased toward runaways due to the magnitude cut preferentially selecting stars with low extinction. Gies & Bolton (1986) observed a population of 36 OB stars selected based on previous identification as high-radial-velocity sources complete to V =8 mag and as faint as V =11 mag. After applying a statistical correction to their observations to account for the sampling bias, they determined a runaway fraction of 2% for early-B and 10% for O stars. Stone (1991) fit a bimodal distribution to the observed velocities of O stars to find raw runaway fractions of 17. 4%, 13.4%, 9.5% for 30, 35, and 40 km s −1 thresholds, respectively. After applying a statistical correction for incompleteness the corrected O star runaway fraction was 46% (Stone 1991). de Wit et al. (2005) examined the O stars identified as field objects by Mason et al. (1998) and found 22 of the 43 (51%) to be runaways based upon their 3D peculiar motions. It should be noted that field populations would be expected to have higher velocities due to their separation from stellar associations and clusters. Lamb et al. (2016) identified a sample of 374 "field" OB stars in the Small Magellanic Cloud. Deriving peculiar radial velocities based upon the radial velocities of nearby gas clouds, they found a lower limit field OB star runaway fraction of 11.3±2.2%, though as noted above, they estimate the true runaway rate may be as much as twice their reported rate if proper motions were included. Based upon a statistical analysis of Gaia DR1 (Gaia Collaboration et al. 2016a) proper motions of Galactic O-Star Catalog stars (version 4.1), Maíz Apellániz et al. (2018) found a runaway fraction of 6.7±1.1%. Seemingly, the only consistent result between these studies is a disjunction between the runaway fractions of O and B stars, with O stars being much higher. These disparate results beg for reexamination with the improved proper motions of Gaia EDR3 data and radial velocities.
Despite the large dispersion in observed runaway fraction, these observational studies conflict with simulations, which predict very low fractions. Two mechanisms for generating runaway stars have been proposed: the binary supernova scenario (BSS), where stars with close binary companions are ejected when the more massive component under goes an asymmetric core-collapse supernova, and the dynamical ejection scenario (DES), where stars are ejected through N-body interactions. Hybrids models have also been proposed to explain the isolated systems of de Wit et al. (2005) through the dynamical ejection of a binary system which then undergoes the BSS. Another hybrid idea involves a supernova in a triple system that ejects a high-velocity close binary (Gao et al. 2019). Studies of the DES mechanism generally find higher proportions of stars achieving runaway velocities and higher fractions of multiplicity among those runaways than the BSS mechanism.
Investigations of the efficacy of the BSS generally rely upon population synthesis models, simulating the evolution of (typically isolated) massive binaries. Such simulations appear highly dependent on their adopted initial model parameters, e.g., initial binary fractions, mass ratios, orbital periods, and whether the simulation code includes mechanisms such as companion mass exchange or N-body dynamical interactions. Assuming an initial binary fraction of 80%, De Donder et al. (1997) found that 5-8% of O stars have runaway velocities >30 km s −1 , and fewer than one-third of those runaway systems are binaries. Modeling only binary systems, Eldridge et al. (2011) found a runaway fraction of approximately 5% for O stars and slightly lower for B stars. Renzo et al. (2019) and Evans et al. (2020) modeled the evolution of isolated binary systems and found BSS-induced O star runaway fractions of 0.5 +1.0 −0.4 % and 0.5-2.0%, respectively. It should be noted that De Donder et al. (1997), Renzo et al. (2019), andEvans et al. (2020) do not consider the additional effect of dynamical ejections in their studies, focusing instead on isolated systems which likely contributes to their lower runaway fractions. Eldridge et al. (2011) acknowledged the likelihood that the DES contributes to the runaway population but predicted only modest increases in the runaway fraction over their strictly BSS predictions.
N-body simulations of stellar clusters predict that dynamically ejected runaways comprise a few percent of the general runaway O star population, similar to the contribution of the BSS. Poveda et al. (1967) modeled 5-and 6-body star clusters analytically and found that 2-17% of stars were ejected with velocities in excess of 35 km s −1 . Using N-body simulations of 10 3.5 M and 10 4.0 M stellar clusters, Oh et al. (2015) found about 30% of the ejected O stars achieve runaway velocities, producing a runaway fraction of 5% along all O stars. One of the few studies that approaches the observed rate of runaways is the DES simulations of Perets &Šubr (2012) that predicted O runaway fractions of 5-20% and B-star runaway fractions of one-third to one-sixth that of the O star fraction.
It appears likely that both the BSS and DES contribute to the population of runaways. Indeed, some observational studies have identified runaway star candidates from clusters that are too young to have experienced their first supernova (e.g., Gvaramadze & Bomans 2008;Roman-Lopes 2013;Drew et al. 2018) making DES the only likely mechanism for producing these systems. Similarly, studies have identified high-proper-motion star+pulsar pairs with anti-parallel velocity vectors that suggest a common co-located origin, making them probable BSS-generated runaways (Hoogerwerf et al. 2000(Hoogerwerf et al. , 2001. Given that each mechanism generally produces runaway fractions of a few to several percent in simulations, combinations of the two mechanisms can only reproduce the low end of the observed runaway fractions, i.e., <10%.

Goals of this Investigation
Here and Paper I we present the results of complementary studies to Peri et al. (2012), Peri et al. (2015), and Bodensteiner et al. (2018). Whereas those works began with catalogs of runaways and bright OBA stars, respectively, and searched for SBNe, we begin with the comprehensive K16 and supplementary Jayasinghe et al. (2019) catalogs of SBNe to determine their spectral types (Paper I) and kinematic properties (this work). With the high-quality measurements from the Gaia Early Data Release 3 (EDR3) (Gaia Collaboration et al. 2016b, it is possible to locate individual stars within the Galactic Plane to distances of nearly 8 kpc, greatly extending the accessible volume over previous stellar catalogs such as Hipparcos (van Leeuwen 2007) or Tycho-2 Høg et al. (2000). Previous measurements of the velocities of SBNe (van Buren et al. 1995;Bodensteiner et al. 2018) have measured proper motions without accounting for the Galactic rotation curve and often without correcting for solar peculiar motion (e.g., K16). While this is reasonable approximation when within the solar neighborhood, at distances of >few kpc these corrections can be on the order of 10-30 km s −1 , similar to the runaway velocity threshold. In this work, we account for both the solar peculiar motion and the expected contribution of the Galactic rotation curve to calculate peculiar velocities within each star's local standard of rest.
The majority of stars in the K16 catalog lack radial velocity measurements, precluding the calculation of true 3D peculiar motions. However, simulations of SBNe suggest that their arcuate shapes become increasingly reniform (circular) as the inclination angle increases beyond ≈45 • from transverse (Acreman et al. 2016;Meyer et al. 2016). As SBNe in the K16 catalog were identified on the basis of an arcuate morphology, this suggests that, if the star-ISM differential is created by the stellar peculiar velocity, the largest components of the peculiar velocity will be transverse and not radial. For this reason we adopt and defend the proposition that the two-dimensional peculiar velocity (v 2D ) calculated from proper motions is a reasonable approximation of the full three-dimensional velocity (v 3D ) for most SBN central stars.
In this work, we calculate v 2D for a sample of SBN central stars and a comparison sample of O stars from Galactic O-Star Catalog (Maíz Apellániz et al. 2013 to determine their velocity distribution and runaway fractions. As a comparison sample, the GOSC contain massive stars earlier than B0 selected for completeness without regard to other properties, whereas the SBN sample is selected on the basis of infrared nebular morphology. Furthermore, the GOSC contains exclusively O stars while the SBN sample contains both O and early B stars-predominantly the latter. If the kinematic properties of O and early B stars differ systematically (e.g., as suggested by Kobulnicky et al. (2016) and Banyard et al. (2022) with regard to O stars having a higher multiplicity fraction than B stars), then the comparisons between the SBN and GOSC samples may be regarded as approximate but imperfect.
Our goals include understanding the origins of runaway stars, the physical mechanisms responsible for production of SBNe, and the role of multiplicity among SBN stars and O stars. If the SBN central stars show an excess of runaways relative to a general O star sample, this would support the idea that stellar peculiar velocities are responsible for the star-ISM velocity differential. However, if the peculiar velocities are the same as the general O-star population, this may indicate that either "in situ" shocks comprise a larger portion of the SBNe population or that ISM density is a more important factor in the formation of an SBNe than the velocity differential. If runaway stars are predominantly single relative to the field O population then the BSS mechanism is the favored one. However, if runaways are commonly multiple-star systems, then the DES mechanism must be responsible for their large velocities.
Section 2 describes the selection of sample sources, the methods and assumptions of the kinematic computations, and spectroscopic reduction techniques for GOSC runaway O stars and SBN central stars observed. Section 3 presents results on the GOSC O-star runaway comparison sample, including the distribution of v 2D , runaway fraction, and multiplicity fraction for a subsample of 455 GOSC sources. Section 4 presents two-dimensional and three-dimensional velocities of the SBN central stars sample, runaway fractions, and the distribution of alignments between the stars' velocities and the orientation of their nebulae. In Section 5 we conclude with a discussion of these results and implications for production of SBNe and runaway stars.

Source Selection
The stellar sources for this work were selected from the K16 catalog of 709 morphologically identified SBNe candidates and their central stars, supplemented by an additional 1 187 citizen-science-identified stars from Jayasinghe et al. (2019). We extracted Gaia EDR3 astrometric quantities for these 896 sources as matched by position using a 1.5 search radius-twice the angular resolution of the mid-IR images that comprise the search dataset-from the position of the mid-infrared point source identified as the probable central star of each nebula. Positional matching resulted in 797 bowshock central stars having a viable EDR3 counterpart. Those without a match are likely objects where the central star is too faint at optical wavelengths to appear in the EDR3. In 19 instances positional matching resulted in two candidate central stars. For the majority of such cases the ambiguity was resolved by selecting the source with the optical magnitude most consistent with the infrared-identified point source IR magnitudes. Of the stars having an EDR3 entry, only 582 have a physically meaningful positive parallax. We removed 78 stars having a renormalized unit weight error parameter (Lindegren et al. 2018, RUWE) >1.4, as these have astrometric solutions showing largerthan-expected residuals-possibly as a consequence of binarity-resulting in less reliable proper motions. Finally, we retained only the stars having parallax and proper motion measurement:error ratios greater than 4:1, leaving 267 objects having well-measured astrometric parameters. The 16th/50th/84th percentiles of the Gaia G-band magnitudes of selected stars are 9.5/12.6/14.8, so the majority of the sources are considered bright targets.
We created a field O star (non-bowshock) comparison sample from the 590-source GOSC, as compiled in Maíz Apellániz et al. (2016) using the same selection criteria as for the SBN central stars, resulting in 455 O stars comprising the GOSC comparison sample. Figure 1 plots the Galactic locations of the 455 GOSC O stars (blue x's) and 267 SBN central stars (red dots) in each sample using a polar representation with the Sun at center and Galactic Center toward right. The vast majority of the SBN central stars lie in the inner Galaxy, with a comparable number in the first quadrant (112 objects) and fourth quadrant (117 objects). Nearly all SBN objects lie within 1 • of the Galactic Plane. The mean distance is larger for SBN central stars (3.3 kpc) than for GOSC O stars (2.3 kpc), a likely consequence of the infrared selection criteria for the former and optical selection criteria for the latter.

Kinematic Calculations
We computed the Galactic position and peculiar velocity (v 2D from the two proper motion components) for each source relative to its local standard of rest based upon the EDR3 proper motions and parallax-based distances. Uncertainties on parallaxes and proper motions are, on average, factors of 3-5 smaller in the EDR3 compared to the earlier DR2 dataset. We adopted inverse parallax values for the distances. The calculation of velocities is based upon the matrix equations of Johnson & Soderblom (1987) in conjunction with an adopted solar motion relative to the Local Standard of Rest (LSR) and an adopted Galactic rotation curve. For a complete description of this technique and the associated matrix transformations, see Appendix A of Randall et al. (2015). To transform from equatorial to Galactic coordinates, we adopted the Reid & Brunthaler (2004) Galactic Center coordinates (α GC = 17 h :45 m :37. s 224s, δ GC = −28 • :56 :10. 23) and Galactic North Pole coordinates (α GNP = 12 h :51 m :26. s 282, δ GN P = 27 • :07 :42. 01). We adopted a Solar Galactocentric distance of 8.15 kpc (Reid et al. 2019) and the Irrgang et al. (2013) Model II Galactic rotation curve. The choice of rotation curve makes very little difference to the calculated velocities because all modern rotation curves are nearly flat within a few kpc of the Sun where most of our sample stars lie. The calculation of peculiar stellar velocities, however, is somewhat sensitive to the adopted Solar motion relative to the LSR. We tried both the solar peculiar motions of Schönrich et al. (2010)  We estimated the uncertainties on the peculiar velocities using a 1000-iteration Monte Carlo analysis, taking into account the EDR3 uncertainties on parallax and the two proper motion components, along with the uncertainties on the solar motion relative to the local standard of rest. This method produces distributions of peculiar velocities in U, V, and W. We then added in quadrature the uncertainties of each peculiar velocity distribution to determine the uncertainty on v 2D . We did not consider systematic uncertainties on the EDR3 motions, as these are estimated to be 0.004 mas yr −1 , negligible compared to the random uncertainties and small compared to the earlier Gaia DR2 systematic uncertainty of ∼0.010 mas yr −1 (Gaia Collaboration et al. 2021). Likewise, the systematic error on the EDR3 parallaxes is estimated at ∼0.010 mas, negligible for the purposes of this work (Lindegren et al. 2021). Uncertainties on EDR3 parallaxes may be up to a factor of two larger than the published values as a complex function of location on the sky, color, and magnitude (Lindegren et al. 2021;Maíz Apellániz et al. 2021;Cantat-Gaudin & Brandt 2021), but we do not consider these owing to the evolving understanding of Gaia systematics as of this writing. Proper motions of bright (G<13 mag) stars may also be systematically rotated with respected to fainter targets (Cantat-Gaudin & Brandt 2021), so we performed an analysis with and without their recommended corrections using their code and concluded that these are insignificant to the conclusions of this work.
For those sources with two or more radial velocity measurements from our spectroscopic survey, we calculated a 3-dimensional peculiar velocity (v 3D ) using the weighted mean radial velocity measurement and adopting the RMS of those measurements as the associated uncertainty. This mean is used to minimize the effects of potential radial velocity deviations from binary (or higher-order multiple) systems, both identified and unidentified.
Adding signed quantities in quadrature (i.e., the two components of proper motion) always contributes positively to the final value. We correct for this positive bias by subtracting in quadrature the velocity uncertainties (σ i ) from the calculated velocities (v i ) for each star, For the selected sources with small fractional uncertainties, this bias correction is small relative to the uncorrected velocities. The median positive bias correction is 1.7 km s −1 for the GOSC stars and 1.8 km s −1 for the SBN stars. In order to test and validate our code, we replicated the Sperauskas et al. (2016) analysis of the peculiar motions of 1088 K-M stars in the solar neighborhood using their set of assumptions and Galactic rotation curve. The average absolute difference between their calculated peculiar velocities and ours is 1.1±1.0 km s −1 .

Spectroscopic Observations of GOSC O Stars
In addition to the optical spectra acquired previously on 84 SBN stars for Paper I, we obtained spectra of 25 stars from the GOSC, selected on the basis of having large v 2D (in excess of 25 km s −1 ), high source brightness (broad range of 7.2 V 11.2 mag), and positive declinations suitable for observation from the northern hemisphere. No sharp magnitude cut was applied so that the sample is not biased toward more luminous (hence, binary) stars. These sources also met the astrometric selection criteria described above. We obtained 119 new optical spectra using the Wyoming Infrared Observatory (WIRO) 2.3 meter and Apache Point Observatory (APO) 3.5 meter telescopes. Each star in this sample was observed on at least three separate nights, with thirteen stars observed more than five nights. The interval between observations varied between one and 84 days with a mean of 30 days and a median of 32 days. The mean time baseline from first to last observation was 107 days with a median of 118 days. All but one star had a minimum baseline of 57 days with 17 of the 25 stars having a baseline longer than 90 days.
We acquired 81 spectra on the nights of 2019 July 29, August 16, 28, and 19, and September 4 using the WIRO LongSlit spectrograph. Observational methodologies for WIRO spectra closely follow those described in Paper I. During each observation, sources were positioned within a 1. 2×120 slit oriented north-south. A 2000 line mm −1 grating yielded a spectral resolution of 1.26Å FWHM (R≈4000) with a spectral coverage of 5400-6800Å. Exposure times ranged from 1×30 s to 2×300 s, depending on source brightness and seeing. CuAr arc lamp exposures were taken before and after each exposure, yielding wavelength solutions with a typical RMS of 0.015Å, or 0.8 km s −1 at 5800Å. The wavelength calibration is estimated to be precise to 6 km s −1 .
We acquired 38 spectra using the Double Imaging Spectrograph (DIS) on the APO 3.5 meter telescope over seven nights: 2019 May 24, 2019 July 3 and 9, 2019 September 4 and 25, and 2019 October 22 and 29. The 1200 line mm −1 grating of the red spectrograph arm yielded a reciprocal dispersion of 0.58Å pix −1 over the 5700-6900Å wavelength range. Depending on source brightness, exposure times ranged from 1×20 s to 2×120 s. The instrument HeNeAr arc lamp supplied wavelength calibration to an RMS of 0.06Å. Instrument rotation during the night produces wavelength shifts of up to 0.3 pix, which were removed using periodic HeNeAr arc lamp exposure such that the wavelength solutions are estimated to be precise to about 5 km s −1 .
Spectra were extracted and reduced using standard longslit techniques within IRAF (Tody 1986). Flat-fielding was performed using quartz lamp dome flats for WIRO spectra and internal quartz lamp flats for APO spectra. Where multiple spectra were taken in a single target, nightly spectra were combined into a single master spectrum yielding a final signal-to-noise ratio (SNR) of at least 65 for APO spectra and 80 for WIRO spectra in the vicinity of the stellar He I 5786Å line used to measure radial velocities. Final spectra were continuum normalized and Doppler corrected to the Heliocentric velocity frame. Radial velocities were extracted using the methodology described in Paper I. Briefly, the Heliocentric velocity for each spectrum was re-zeropointed by shifting it (typically 3-15 km s −1 ) to place the centroids of the interstellar Na I D λλ5890, 5896 lines at the mean velocity for the ensemble of each star. Using a fitting code based around the IDL routine MPFIT Markwardt (2009), we measured radial velocities from each spectrum by fitting Gaussian functions to the He I 5786Å line. This fitting code fixed the Gaussian width and depth to the mean from the ensemble after manually rejecting outliers and solved for the best line center and uncertainty. This technique allowed us to measure robust relative radial velocities but, since our data lack radial velocity standards, these should not be treated as absolute radial velocities due to potential systematic offsets relative to the Heliocentric frame. We estimate any such absolute systematic effects to be small, <8 km s −1 , based on calibrations of other surveys using the same instrument (Kobulnicky et al. 2014).

Runaway Fraction of the General O Star Population
Figure 2 presents a histogram of the distribution of v 2D for the 455 GOSC stars, using a bin width of 2.5 km s −1 . Vertical dashed lines mark the two fiducial thresholds for runaways at v 2D >25 km s −1 and v 2D >30 km s −1 , respectively. Five systems lie off the right edge of the plot at v 2D >100 km s −1 . The distribution in Figure 2 is similar to the peculiar tangential velocities of Tetzlaff et al. (2011, Figure 4) having a peak near 8 km s −1 with a long tail extending out to 100 km s −1 . Tetzlaff et al. (2011) do not apply a positive bias correction so their distributions may be biased to slightly larger velocities. The median/mean velocity uncertainties are 1.8/2.6 km s −1 .  Of the 455 GOSC sources, 81 +14 −8 , have v 2D >30 km s −1 , yielding a runaway fraction of 16 +3 −1 %. For this work, we adopt a v 2D runaway threshold of 25 km s −1 as our primary runaway star criterion, for two reasons. Tetzlaff et al. (2011) analyzed 7663 stars within 3 kpc of the Sun and found the peculiar tangential velocities (our v 2D ) could be described by a two-component Gaussian distribution, with low-and high-velocity components peaking at 7 and 22 km s −1 , respectively. These low-and high-velocity distributions intersect at 20 km s −1 , with the low-velocity component approaching zero amplitude at ≥25 km s −1 . Therefore, adopting 25 km s −1 as our v 2D threshold would select only members of this high-velocity group. Furthermore, a v 3D threshold of 30 km s −1 corresponds to a v 2D threshold 25 km s −1 ( 3/2≈30/25). After applying the positive bias correction to the data, our v 2D threshold of 25 km s −1 yields 102 +14 −15 of 455 GOSC O stars as runaways (22 +3 −3 %). Our GOSC runaway fraction is notably higher than that found by Maíz Apellániz et al. (2018) who analyzed 594 GOSC stars (including some previously unpublished stars) using Gaia DR1 proper motions. However, Maíz Apellániz et al. (2018) used a statistical analysis based on 2D proper motions corrected for solar motion, rather than a true peculiar velocity. This method was used for computational simplicity and to avoid the use of DR1 parallaxes which did not include bright sources and had large relative uncertainties. The statistical cut for their runaway candidates was empirically determined using the results of Tetzlaff et al. (2011). Maíz Apellániz et al. (2018 found 49 runaway candidate O stars out of 594 sources, yielding a runaway fraction of 8.2%, considerably lower than our 22±3%, as well as lower than historical literature values for O stars. Opting to use the v 2D >30 km s −1 fraction of 17 +4 −1 % reduces but does not resolve the inconsistency. One potential reason for their lower runaway rate may be their analytical method being insensitive to stellar distance. In our analysis 77 (17%) of the GOSC sources have distances greater than 3 kpc. At greater distances large transverse velocities have proportionally smaller proper motions. Because of their focus on proper motion outliers, their analysis may exclude high-peculiar-velocity stars at large distances. Of the Maíz Apellániz et al. (2018) runaway candidates our analysis recovered all 30 that met our selection criteria. Our new field O star runaway fraction of 22±3% is consistent with most previous observations of O stars as outlined in Section 1, including the raw fraction found by Stone (1991) and the Blaauw (1961) O stars, but substantially lower than the Stone (1991) corrected fractions of 46% and lower than the 55% reported by de Wit et al. (2005). Table 1 lists the computed stellar kinematic properties for each of the 25 spectroscopically studied GOSC runaway candidates. Columns 1 and 2 list the GOSC catalog identification number and a common alias, respectively. Column 3 is the EDR3 RUWE. Columns 4 and 5 are the calculated v 2D and 1σ uncertainty. Column 6 is the mean Heliocentric radial velocity of the 3-6 individual measurements, column 7 is the average velocity uncertainty, and column 8 is the RMS of the 3-6 radial velocity measurements. Column 9 gives probability P (χ 2 , ν) for each star's set of radial velocities being consistent with the null hypothesis (no velocity variability) using a chi-squared test with ν degrees of freedom, where ν is the number of radial velocity measurements minus one. Sources with P (χ 2 , ν)<0.05 are inconsistent with the null hypothesis of a single-star system and are interpreted as candidate multiple-star systems. Column 10 lists the spectral type and luminosity class from the literature, and column 10 gives the dates of observation for each star in a modified Heliocentric Julian Date format (HJD−2,458,000). Figure 3 displays a histogram comparing the distributions of P (χ 2 , ν) for the 25 spectroscopically studied GOSC O stars and the 72 SBN central stars having multi-epoch radial velocity data presented in Paper I. The distributions of P(χ 2 , ν) are similar for the two populations, showing a roughly flat distribution as expected for a population of single-star systems but with a sharp peak at P (χ 2 , ν)<0.05 indicating a population of candidate multiple-star systems. Of the 25 GOSC systems, 14 (56±16%) have P (χ 2 , ν)<0.05 compared to 27 of the 72 SBN central stars (38±6%). If we were to consider only the 52 SBN stars meeting our astrometric selection criteria, then 16 (31±7%) have P(χ 2 , ν)<0.05. Assuming the average occupancy of those bins with P(χ 2 , ν)>0.05 to be an estimate of the false positive rate for each sample yields corrected multiplicity numbers (fractions) of 13.4 GOSC stars (53±14%) and 25.6 SBN central stars (36±6%). If we consider only the O stars among the SBN central stars, 3 of 8 (38±22%) show evidence of multiplicity. These results show that both the full OB-star SBN sample and O-star SBN sample are consistent with the runaway GOSC O-star multiplicity fraction within 1σ. The apparent multiplicity fraction of the GOSC and SBN samples are formally consistent, despite the fact that early B stars dominate the SBN sample , and B stars have been shown to exhibit a lower multiplicity fraction than O stars (Kobulnicky et al. 2014;Banyard et al. 2022). This bias should augment the disparity in spectral type between the samples. The consonance between SBN subsamples suggests that SBN central stars, despite consisting mostly of early B stars, may have larger multiplicity fractions than GOSC O stars. Or, this result could be interpreted to mean that binary systems are more effective in generating observable stellar bow shock nebulae, since the SBN sample is constructed exclusively by the presence of these nebular infrared features.

Multiplicity Fraction of Runaway GOSC Stars
We considered the possibility that quasi-periodic line profile variations attributed to non-radial pulsational modes and/or stochastic atmospheric fluctuations among early type stars (rather than multiplicity) produce some of the observed variability and artificially inflate the multiplicity fraction of the GOSC sample. Fullerton et al. (1996) found that 23 of 30 O stars-and all of the supergiants in their sample-exhibited significant line profile variations in resolution R≈20,000-30,000 optical spectra. Some of these would be large enough to mimic the signature of binarity seen in our R=4000 spectra. However, in a spectroscopic survey of 128 O and early B stars that yielded orbital solutions for 48 systems using at least 12 observations per system, only six stars were found to exhibit aperiodic line variations attributed to atmospheric effects (Kobulnicky et al. 2014). All six were supergiants. Among the present sample, 7 of the 14 GOSC stars showing velocity variability are supergiants. If we conservatively suspect all of the supergiants of being dominated by pulsational variations and remove them, this leaves 7 of 25 (28%) as multiple-star candidates. A more modest correction based on the Kobulnicky et al. (2014) survey suggests removing 2-3 velocityvariable supergiants as potential contaminants. These considerations lead us to adopt a revised multiplicity fraction of >28-48% for GOSC runaway O stars.
All of these fractions should be considered minimum multiplicity fractions considering the short observational time baseline for the GOSC spectroscopic observations. The data on the SBN central stars presented in Paper I, with its longer observational time baselines, was able to identify all five previously recognized spectroscopic binary sources. Among the GOSC spectroscopic sample, five stars (GOSC IDs 27,206,234,340,365) are identified in the literature as spectroscopic binaries, primarily by Pourbaix et al. (2004). Our spectroscopic program identified three of these (27, 206, and 365) as having P (χ 2 , ν)<0.05 but did not identify the two other literature binaries. GOSC 340 (HD 192281) was originally identified as a binary system by Barannikov (1993) with a period of 5.48 days and a semi-amplitude of 16.8±2.4 km s −1 (though Maíz Apellániz et al. (2019) note this has not been subsequently confirmed). The average velocity uncertainty for this source in our dataset is 13 km s −1 , of a similar magnitude as the semi-amplitude. The orbital period of GOSC 234 (HD 15137) is 28.6±0.4 days McSwain et al. (2007). The irregular observational cadence of our study was 56.0, 1.0, and 60.7 days between observations, approximately two periods. This unfortunate alignment of cadence and orbital phase explains why our data do not reveal its multiplicity. Including these two literature binaries would raise the raw (uncorrected for line profile variability contamination) multiplicity fraction of our GOSC sample to 16 out of 25 sources-64%-comparable to multiplicity fractions of other OB star samples (Kobulnicky et al. 2014;Banyard et al. 2022). Our data reinforce the conclusion that multiple-star systems are common, not just among OB stars in general, but among high-velocity runaways as well.

Kinematics of the Central Stars of Stellar Bowshock Nebulae
We divided the SBN central stars that met our astrometric selection criteria into three subsamples. The first group is the "Multiple-Star Candidates": 16 sources, selected on the basis of having P (χ 2 , ν)<0.05. The second group is the "Single Star Candidates": 36 sources having P (χ 2 , ν)≥0.05. There are likely to be multiple-star systems in this group that were not identified in our radial velocity survey. However, the fraction of binaries in this sample will be lower than a typical SBN central star sample. The third group, termed "Group 3", consists of the 215 sources that met our astrometric selection criteria but do not have measured radial velocities. This sample is likely the most representative of SBN sources in general based upon the multiplicity cuts for the first two groups. Table 2 provides identifications and tabulated astrometric data for the 267 SBN central stars meeting our astrometric selection criteria. Column 1 is either the K16 catalog identification number from Kobulnicky et al. (2016, numbers 1-709) or an identification index that includes the 187 additional stars drawn from Jayasinghe et al. (2019, numbers 710-896). Column 2 is a common alias for the star. Columns 3 and 4 are the right ascension and declination of the star, in degrees. Column 5 is the EDR3 ID number. Column 6 gives the Gaia G magnitude. Column 7 is the environment class from K16 specifying the local environment of the nebula as isolated (I; 155 instances), facing an H II region (FH; 50 instances), within an H II region (H, 30 instances), or facing a bright-rimmed cloud (FB, 30 instances). An FB designation often entails an H designation, so some objects have a compound environment class (e.g., FB/H). Columns 8 and 9 give the EDR3 parallax and uncertainty, in milliarcsec. Columns 10-13 give the EDR3 proper motion in right ascension and declination, in milliarcsec yr −1 , and their associated uncertainties. Column 14 is a code designating each star as part of the Multiple Star Candidates group (1), Single Star Candidates group (2), or the Group 3 (3). Column 15 is the spectral type, either from Paper I where sources were classified into three broad categories O/OB/B, or from the literature, in which case the spectral type and luminosity class is enclosed in parentheses. Column 16 lists the literature reference for the spectral classifications in parentheses.
The literature spectral types of Table 2 reinforce the conclusions of Paper I that the central stars of SBNe are predominantly OB spectral types. Of the 141 sources with spectral types, 136 are OB stars, three are M stars (possibly the evolved descendants of massive stars), one is a carbon star, and one is a K star. As part of this work, we examined the Wide Field Infrared Explorer (WISE) and Spitzer Space Telescope archival images to re-evaluate local environments of the SBN central stars using the criteria of K16. We cannot determine that the arcuate ISM structures are co-distant with the SBN central stars, only that these structures are co-located on the sky. For the majority of sources (65%), we affirmed the environment class of K16. Of the sources for which we reassign the local environment class, 13 of 56 resulted from changing identifications within complicated regions (e.g., updating a source from an H to an FB where the rim of an interstellar bubble was prominent at 8.0 µm). Nearly half, 24 of 56, involved updating sources from I to FB or FH. These sources were largely 5-15 arcmin from the associated H II region or bright-rimmed cloud, many without sharp boundaries. The remaining 19 sources we reassigned the environment class from I to H, where these stars and H II regions appear at the edge of WISE ATLAS tiles. For the 50 sources in Table 2 drawn from the extended bowshock candidate list of Jayasinghe et al. (2019), the listed environmental classes are new to this work. Table 3 lists kinematic measurements for the SBN central stars. Column 1 is the identification number, as in Table 2. Columns 2 and 3 give the calculated v 2D and associated uncertainty in km s −1 . Where radial velocities were measured for the Multiple Star Candidates and Single Star Candidates, Columns 4 and 5 provide the calculated v 3D and associated uncertainty in km s −1 . The largest source of uncertainty in the v 2D calculations is the distance uncertainty arising from the linear dependence of each velocity component on distance. The uncertainties on the angular proper motions, σ µα and σ µ δ , average <1%, whereas the mean parallax uncertainty is 8%, thereby constituting the largest source of error on v 2D . Uncertainties on v 3D will always be larger than the uncertainties of v 2D due to the additional uncertainty term from the radial velocity measurement. For the multiple-star systems, this additional uncertainty may be significant, ≥6 km s −1 as an observational minimum to as much as 50 km s −1 for some binary systems where the systemic radial velocity is poorly constrained from just a few measurements. Column 6 lists the morphological position angle (P A m ) from the central star to the apex of the SBN, in degrees from North toward East in Equatorial coordinates, as tabulated in K16. This position angle was measured by eye based upon Spitzer Space Telescope 24 µm and WISE 22 µm archival images and carries an estimated uncertainty of 8 • . Columns 7 & 8 are the kinematic position angle (P A k ) of the star's v 2D vector and its associated uncertainty in degrees. Columns 9 & 10 give the difference between the morphological position angle and the kinematic position angle in degrees, ∆ PA =P A m − P A k , and associated uncertainty. Column 11 provides the code, as in Table 2, designating each source as part of the Multiple Star Candidate group, the Single Star Candidate group, or the Group 3 sources which lack radial velocity measurements. Figure 4 plots the calculated v 3D versus v 2D for the 25 systems in the GOSC runaway sample (blue symbols) and the 52 SBN central stars having radial velocity data (red symbols). Circles denote multi-star candidates, while triangles denote single-star candidates. A solid black diagonal line marks the 1:1 relation where v 3D =v 2D (i.e., where a source would have zero peculiar radial velocity). A dashed line marks the relationship for isotropic velocities where the :v 2D ratio of 2:1. The majority of SBN stars, 31 of 52 (60%), lie between the solid and dotted lines and are consistent within one sigma of the dashed line, indicating that their velocities are either entirely within the plane of the sky or consistent with an isotropic velocity distribution. This agrees with the expectation that morphologically selected SBNe central stars have velocities primarily in the plane of the sky. If a star's velocity were primarily radial through a stationary ISM, the morphology of the nebula would appear more circular than arcuate (Meyer et al. 2016). For a minority of sources (40%), particularly in the low-v 2D regime ( 10 km s −1 ), the ratio v 3D :v 2D is greater than 2. These are likely to be either unidentified binary systems or multiple-star systems with a poorly determined systemic velocity, a consequence of having only a few radial velocity measurements. One system, the SBN multiple-star candidate 346 with v 3D =220 km s −1 , lies off the top left of the plot. Upon closer examination, SBN 346 appears less arcuate and less symmetric than the other sources in the K16 catalog. Given its large peculiar radial velocity of 196±9 km s −1 and dubious nebular morphology, we reclassify this as a doubtful SBN candidate. It may be a radial velocity runaway or a multiple-star system with a poorly determined radial velocity.   . Three-dimensional peculiar velocity versus two-dimensional peculiar velocity for multiple star candidates (circles) and single-star candidates (triangles) from the SBN sample (red) and the GOSC comparison sample (blue). Only stars having multiple radial velocity measurements and meeting the astrometric selection criteria are shown. A solid line marks the 1:1 correspondence, while a dashed line marks the 3/2:1 relation expected for isotropic velocities in three dimensions.
For some SBN central stars the calculated v 2D may not be representative of the true star-ISM velocity differential if there exists a significant local flow of interstellar material. Such a local flow will be a significant factor in determining the nebular morphology and the relative star-ISM velocity, most dramatically so if the star's peculiar velocity is small. Winds from H II regions may impart local flow velocities as large as 30 km s −1 (Tenorio-Tagle 1979; Bodenheimer et al. 1979;Castaneda 1988). If a local flow is transverse, the star-ISM velocity differential would also be transverse creating the arcuate morphology without need for a stellar peculiar velocity. This may help explain why stars with peculiar velocities seemingly dominated by radial velocity components exhibit arcuate SBN. This class of objects would include the "in-situ" bowshock candidates where an arcuate nebula points toward a nearby H II region (e.g., our FH environment class). Of the 267 SBN objects in this sample, 50 (19%) have FH environmental classifications, meaning that a non-negligible fraction could be influenced by local ISM flows.

The Runaway Fraction of SBN Central Stars
Figure 5 presents histograms and kernel density estimates (KDEs) for the three SBN central star groups: Group 3 systems (cyan), the Single-Star Candidates (magenta), and the Multiple-Star Candidates (green). Colored curves denote the KDE for each group and the full SBN sample (black dashed curve). An unfilled black histogram and black solid curve designate the distribution and KDE of the 455 GOSC stars, renormalized to facilitate comparison. The KDEs for each population show a distinct peak centered near 9 km s −1 and a long tail out to 80 km s −1 and beyond. Although the KDEs suggest the Single Star Candidates may contain a high fraction of runaway stars, this may be a result of the small number of stars in that group. All SBN groups are consistent with having been drawn from the same parent population with the p-values (from K-S tests) >0.30 for all sample comparisons. However the SBN and the GOSC distributions are different at the p=0.003 level. This difference is mostly driven by the B stars within the SBN sample, as a comparison of just the SBN O stars with the GOSC sample results in p=0.042, marginally consistent with being drawn from the same population. The SBN sample is noticeably shifted toward larger velocities and has fewer objects at small velocities. The 16th/50th/84th percentile velocities of the SBN sample are 6.5/14.6/32.6 km s −1 versus 5.6/11.4/32.6 km s −1 for the GOSC sample. This difference could be attributed to the larger distances for the SBN stars (mean distance of 2300 pc) than the GOSC stars (mean distance of 3300 pc), since transverse velocity inferred from proper motion scales linearly with distance.
From  Among the 142 stars with spectral classifications from the literature or our work (Paper I), 43 (30%) are O stars, 61 (43%) are early B stars, with the remainder (26%) being indeterminate OB stellar types on the O9/B0 boundary (other than the six late-type stars F-M). Of the 43 O stars, 12 +2 −1 have v 2D > 25 km s −1 (28 +4 −2 %). This O-star runaway fraction is larger than the 22±3% observed among the GOSC stars and is consistent with the high end of the observed O star runaway fractions (e.g., Blaauw 1961;Stone 1991;Oh et al. 2015). Among the 61 B stars, 16 +4 −1 (26 +7 −2 %) are runaway star candidates, much larger than the 2% B star runaway fractions measured by Blaauw (1961) and Gies & Bolton (1986

Morphological and Kinematic Position Angles
We used the computed ∆ PA parameter-the difference between morphological and kinematic position angles-to assess the degree of alignment between the SBN morphological position angles and the stars' peculiar velocity vectors. Figure 6 presents a histogram of the ∆ PA distribution for the 267 SBN central stars, using bin widths of 10 • . The Figure shows a peak near ∆ PA =0 • , indicating strong alignment between the morphological axes of the SBN and the peculiar motion vectors. This peak is present regardless of the choice of Solar peculiar motions (see Appendix), but is more pronounced by about 50% when the adopted [U , V , W ] values are used. There is a broad tail to either side, with some sources showing anti-alignment (∆ PA =±180 • ). The distribution can be approximated by two hypothetical populations: a highly aligned population centered near 0 • having a Gaussian width of σ=25 • (dashed curve) and a random (non-aligned) population. The solid black curve depicts the sum of these two components. The "highly aligned population" comprises 31% of the sample, while the "random population" comprises 69%. The overall distribution is similar to that presented in K16 (Figure 12) but contains more sources, uses better astrometric data, and corrects for solar motion relative to the Local Standard of Rest.
The two populations identified on the basis of Figure 6 suggest that, in a minority of cases, the stellar peculiar velocities are responsible for shaping the orientation of the nebulae. In a majority of cases, there is considerable dispersion between the star's velocity vector and the morphological axis of the nebulae. This suggests that local environmental factors, such as local ISM flows or the proximity to an H II region or bright-rimmed cloud undergoing photoevaporation-rather than the stellar motion-dominates the formation of the nebulae in most instances. Notably, in about 12% of the sources, |∆ PA |>100 • , indicating anti-alignment between morphological and kinematic position angles. Such anti-alignment could occur in cases where an OB star is ejected form a nearby H II region and an outflow from that H II region impinges upon the star.   Figure 6 owing to the smaller numbers of sources in the latter two groups. We merge the H and FB designations, as stars inside H II regions often face bright rims lining those cavities and may encounter similar photoevaporative flows. The largest group, the isolated stars (I), shows a similar distribution of ∆ PA to the general population in Figure 6. Of the 155 stars in isolated environments, 115 (74%) are highly aligned, having |∆ PA |<45 • . The two remaining groups, H/FB and FH, both have flatter distributions with small, but not statistically significant overdensities near ∆ PA = 0 • . Assessed another way, the Pearson correlation coefficient between P A m and P A k is r=0.37 (p=10 −6 ) for the 155 environment class I stars, confirming the strong correlation. The 62 FB/H environment class stars have r=0.11 and p=0.36, indicating a weak or no significant correlation. The 50 environment class FH stars also show no correlation (r=0.04, p=0.75). The weak or non-existent correlation for FB/H stars is unsurprising, as we would expect significant small-scale flows in the proximity of bright rimmed clouds and bubbles where photoevaporative flows are common and significant relative to the peculiar velocity of the stars. The lack of correlation for FH stars is also unsurprising. We would also expect flows emanating from H II regions-at velocities that may dominate the relative star-ISM differential in cases where the star's velocity is small. Finally, we stress that our environmental classifications are an imprecise characterization of the stars' true local conditions, based only on angular separations from prominent infrared features. A more quantitative analysis comparing SBNe three-dimensional Galactic positions with those of OB stars and H II regions would improve the reliability of the environmental classifications. Such an effort should now be possible but is beyond the scope of this work.  with the largest velocities also have the smallest ∆ PA , falling within the shaded region. Forty-four of the 155 isolated stars (28±4%), eight of the 62 FB/H stars (13±5%), and ten of the 50 FH stars (20±6%) are runaways. Thirty-eight of the 63 runaway stars (80%) show strong alignment (within 45 • ) between their morphological and kinematic position angles. Among stars with projected velocities <25 km s −1 , there is a large dispersion in ∆ PA . Evidently factors other than peculiar velocity determine ∆ PA for these nebulae, such as local ISM flows that are on the order of the ≈10 km s −1 peculiar velocities measured for this population. Other factors may include density inhomogeneities in the vicinities of spiral arms, molecular clouds, or star clusters that create departures from the smooth circular Galactic rotation curve assumed in computing v 2D .
Among the runaways a few stars show dramatic misalignment, having ∆ PA >120 • . Star 303 (labeled in Figure 8) is among the most distant in the sample (d>10 kpc at Galactic longitude =60 • ), giving it a parallax uncertainty that results in a very large uncertainty on v 2D (σ v2D = 68 km s −1 ). Formally, it is consistent with having a negligible (non-runaway) peculiar velocity. Star 865 (v 2D =108 km s −1 ) is also quite distant (> 5 kpc at =340 • ) but has well-measured kinematics. It appears to be confined within an interstellar cloud or bubble irradiated by an external source of illumination, so the tabulated morphological position angle probably corresponds to the direction of this radiation and not the stellar velocity vector. The nebulae preceding Star 873 is compact and appears isolated from external influences; we have no good explanation for the large ∆ PA =−147 • . Star 880 has a well-defined nebula and well-measured kinematics, yet its morphological orientation is nearly 180 • from its projected velocity. Star 624 faces an H II region, has a velocity just above the 25 km s −1 threshold, and shows anti-alignment at ∆ PA =170 • .

CONCLUSIONS
The Gaia EDR3 astrometric data and new multi-epoch optical spectroscopy has allowed us to investigate the kinematic properties of a large sample of morphologically selected stellar bowshock nebulae and their central stars. Quantities of particular interest for addressing the origins of bowshock nebulae and runaway stars are the stellar peculiar velocities (both 2D and 3D), their multiplicity, as evidenced by the dispersion in radial velocity measurements, and the degree of angular alignment between the SBNe morphological axes and stellar kinematic vectors. For comparison we have computed the same quantities for a control sample of Galactic O stars. Table 4 provides a summary of key statistics derived from this study regarding runaway fractions, multiplicity, and alignment of SBN morphological and central star kinematic axes.
Based on the v 2D velocities the runaway fraction for SBN central stars is 24 +9 −7 %, nearly identical to the 22 +3 −3 % runaway fraction for GOSC stars. The SBN runaway fraction is consistent with the high end of simulations (Oh et al. 2015) and observations (Blaauw 1961;Stone 1991), but generally larger than most estimates. The SBN sample consists primarily of early B stars, and so the runaway fractions are much larger than other estimates of B stars which lie near 2% (e.g., Blaauw 1961;Gies & Bolton 1986). We conclude that while the runaway rate of SBN central stars is likely enhanced relative to the general OB population, runaway stars do not comprise a majority fraction of SBN central stars. The median v 2D of our sample is 14.6 km s −1 , far below runaway thresholds and contrary to some common assumptions that bowshock central stars have large peculiar velocities. Both the Multiple-Star Candidates group and the Single-Star Candidate group show similar runaway fractions, indicating that the mechanism responsible for their large velocities apparently acts on binary systems with equal efficiency as single stars. This observation would seem to favor the DES (N-body interactions in clusters) over the BSS (supernova ejection) scenario.
Four of the 13 SBN runaways that have multiple spectroscopic measurements show evidence for binarity, yielding a minimum multiplicity fraction for runaways of 30 +16 −15 %. Among the GOSC runaways in Table 1, this minimum multiplicity fraction is >28% after corrections for false positives arising from potential line profile variations caused by atmospheric pulsation. Both of these fractions are much larger than the ∼few percent proposed by most pure BSS ejection simulations (e.g., Eldridge et al. 2011;Renzo et al. 2019). DES ejection simulations predict 5-20% (e.g., Oh et al. 2015;Perets &Šubr 2012), in better agreement with-but still markedly lower than-the observed fractions for SBN central stars and Galactic O stars. If SBN central stars are representative of the population of field OB stars, as is suggested by the consonance between the samples' multiplicity fractions and runaway fractions, this suggests that the DES contributes more heavily to the field OB population than BSS. The significant numbers of multiple systems ejected from young stellar clusters (Maíz Apellániz et al. 2022) adds evidence to this scenario. Alternatively, given the high fraction of SBN central stars in isolated environments with low peculiar velocities ( 14 km s −1 ), this may provide evidence that a fraction of massive stars may form in relative isolation. However, without attempts to trace these stars back to known star forming regions and OB associations, this result is largely conjecture and begs future follow-up studies that attempt to link isolated massive stars to their birth locations.
About one-third (31%) of SBNe have morphologies that are well-aligned with the motion of their central stars-a "highly aligned population". The majority (69%) show weak or no alignment-"the random population". The degree of alignment increases with peculiar velocity (Figure 8). The relative lack of alignment for low-velocity stars implicates the effects of local ISM flows particularly in the vicinity of energetic H II regions. Indeed, the sources showing low degrees of alignment tend to be those facing or within H II regions or facing bright-rimmed clouds where local flows are expected to play a role (Figure 7). However, some seemingly isolated nebulae also show low levels of alignment. Given the low stellar peculiar velocities (median of 14.6 km s −1 for SBN stars), these flows do not need to be especially fast to influence the orientation of the nebula. The flattened distribution of ∆ PA for weakly aligned sources (Figure 6), may indicate that the majority of SBN central stars are isotropically distributed "walkaways" acting as the "interstellar wind vanes" (Povich et al. 2008). Enhancements in gravitational potential, such as along spiral arms, may also play a role in creating local deviations from the smooth circular Galactic rotation curve assumed in our analysis, at the level of 5-10 km s −1 . Given the consistency of SBN central star properties to those of O stars in the field, stellar bowshock nebulae are a useful infrared indicator of the presence of an OB star, particularly for identifying field populations that lie far from young clusters and OB associations.

ACKNOWLEDGMENTS
An early version of this manuscript appeared as Chapter 3 in the PhD dissertation of W.T. Chick . We thank Jesús Maíz Apellániz and Doug Gies for helpful comments on an early version of this manuscript. We thank students Logan Jensen, Harrison S. Leiendecker, Jacob N. McLane, Evan Haze Nunez, Jason A. Rothenberger, and Ashley N. Piccone for taking observations that contributed to this work. Our team is grateful for support from the National Science Foundation through grant AST-1412845, AST-1411851, AST-2108347, and REU grant AST-1063146, as well as NASA through grant NNX14AR35A. 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.

A. A NEW APPROACH TO COMPUTING LOCAL STANDARD OF REST VELOCITY
Astrometric and radial velocity measurements in the Heliocentric frame of reference require corrections for Solar motion relative to the Local Standard of Rest before conversion to a rotating Galactic reference frame. Many works have calculated [U , V , W ], the solar motion relative to the LSR toward Galactic center, toward the direction of prograde Galactic rotation, and toward Galctic North, respectively. A wide range of values has been reported-U =6-12 km s −1 , V =3-26 km s −1 , W =5-9 km s −1 -even as improved kinematic data, larger stellar samples, and better understandings of systematic biases by stellar age have emerged (e.g., Mayor 1974;Dehnen & Binney 1998;Bobylev & Bajkova 2007;Aumer & Binney 2009;Schönrich et al. 2010;Golubov et al. 2013;Ding et al. 2019). The sample under consideration here, OB stars, represent a young stellar population that likely suffer from some of the same biases (i.e., a low velocity dispersion relative to older populations) as discussed in these works. Our attempts to adopt recent estimates for [U , V , W ] (e.g., Schönrich et al. 2010;Ding et al. 2019) yielded stellar kinematic vectors that showed physically implausible trends by Galactic longitude. This led us to pursue using the SBN sample to derive a new estimate of the Solar motion parameters. Our intent was not a repudiation of other recent measurements of [U , V , W ], which are based on more extensive analyses of larger stellar samples, but rather an ad hoc measurement based on a fundamentally new approach that yields physically plausible results for the SBN OB stars. Our results turn out to be consistent with the range of other recent Solar motion velocities.
Our initial analyses using either the Schönrich et al. (2010)  .26] km s −1 solar velocities produced systematic trends in stellar kinematic position angle with Galactic longitude. Using the former set of Solar velocities yielded SBN central star velocities that were systematically directed toward smaller longitudes in the first quadrant and larger longitudes in the fourth quadranta result that is physically implausible. Under the latter assumption the opposite systematic trend obtained-also implausible. An incorrect Solar motion vector can explain these systematic effects.
The choice of stellar population used to infer [U , V , W ] changes the computed values, most dramatically so with V (e.g., see discussion in Aumer & Binney 2009;Golubov et al. 2013). Young populations (OB stars) have colder kinemtics (smaller velocity dispersion) than older populations which have had time to be heated within the Galactic potential. For the analyses in this work we chose to estimate an ad hoc solar motion using our SBN dataset itself by applying a set plausibility constraints that serve as boundary conditions for our estimation.
1. The kinematic vectors of SBN central stars should not correlate with Galactic longitude. Although our sample lies mostly in the first and fourth quadrants, their peculiar motions should have no preferred orientation by quadrant in order to be consistent with the expected isotropic velocities.
2. The dispersion between the nebular morphological position angles and the central star's kinematic position angles should be minimized. Although one of the purposes of this work is to explore the extent to which morphological and kinematic position angles align, such an alignment is already apparent even without correction for LSR Solar motion. Taken as a whole, the sample distribution of SBN nebular position angles should be maximally consistent with the distribution of central star kinematic position angles. Figure 9 illustrates the distributions of kinematic position angles for SBN central stars resulting from Solar motion parameters [U , V , W ]= [5.5, 7.5, 4.5] km s −1 , adopted as values that best satisfy the above bounday conditions. The upper panel plots Galactic longitude versus kinematic position angle in Galactic coordinates. Here, it is clear that the sample preferentially lies in the first and fourth Galactic quadrants with few objects in the second and third quadrants. There is no probable correlation (p=0.36) between longitude and position angle, as expected for isotropic stellar motions. The lower panel shows histograms of the bowshock nebular morphological position angles in Galactic coordinates (red filled histogram) and the central star kinematic position angles (dashed black histogram). A Kolmogorov-Smirnov (K-S) test shows that the two distributions have a 52% probability of being drawn from the same parent distribution. Both histograms show excesses near Galactic position angles 90 • and 270 • , i.e., along the Plane. A similar excess was noted for the SBN parent sample morphological position angles in Kobulnicky et al. (2016). We understand this excess to be a selection effect rather than a violation of isotropy. High-velocity OB stars having vectors orthogonal to the Plane will traverse the ≈80 pc scale height (Clemens et al. 1988) of the inner disk's molecular layer even in their short lifetimes, making them less likely to generate detectable bowshock nebulae. Coversely, OB stars having velocity vectors parallel to the Plane will remain within this higher density layer longer and produce observable bowshock nebulae for a greater fraction of their lifetimes. Our adopted [U , V , W ] not only provide stellar motions that pass the self-consistency checks outlined above, they also lie within the range of values obtained by other (more traditional) types of analysis, despite being grounded in a fundamentally new methodology. U =5.5 km s −1 lies at the low end of-but consistent with-published ranges (see Table 1 of Ding et al. 2019, 5.5-11.7 km s −1 ). V =7.5 km s −1 falls near the median of published values (3-14 km s −1 , with a couple estimates near 20 km s −1 ). Our W =4.5 km s −1 is considerably smaller than the narrow range of 6.5-7.5 km s −1 generally produced by other studies. Values larger than 5.0 km s −1 produce significantly poorer agreement with our stated boundary conditions. It may be the case that the OB central stars of SBN have a systematically smaller z component relative to other ( The two are consistent with having been drawn from the same parent population. te-(1) GOSC source identification number, (2) common alias, (3) Gaia EDR3 re-normalized unit weight error, (4) projected two-dimensional peculiar velocity v2D, (5) uncertaint on v2D, (6) mean Heliocentric radial velocity, (7) mean radial velocity uncertainty, (8) radial velocity RMS, (9) probability of the null hypothesis; probabilities <0.05 are considere multiple-star candidates, (10) literature spectal type and luminosity class, (11) Dates of observation for radial velocity data in Heliocentric Julian Date minus 2,458,000.                           Multiplicity fraction among 25 runaways >28% 16th/50th/84th percentiles of v2D 5.6/11.4/32.6 km s −1