© ESO 2022

1. Introduction

Dynamical interactions between three or more bodies in stellar clusters and supernova (SN) explosions in binary systems can lead to stars moving at high speed through the Galaxy (Zwicky 1957; Blaauw 1961; Poveda et al. 1967). The ejected stars with velocities higher than 30 km s⁻¹ are called runaway stars (Hoogerwerf et al. 2001) and can be easily found with Gaia astrometry (Maíz Apellániz et al. 2018) while the slower ones, called walkaway stars, are more common but more difficult to identify (Renzo et al. 2019a). Ejections favor high-mass stars over low-mass ones and tend to occur early in the cluster life (Oh & Kroupa 2016). Hence, such events can have a large impact on the massive-star mass function (Pflamm-Altenburg & Kroupa 2006) and, in extreme cases, they can produce “orphan clusters” where the most massive stars have disappeared from the cluster altogether and the present-day mass function (PDMF) is capped with respect to the initial mass function (IMF, Oh et al. 2015; Villafranca II, Maíz Apellániz et al. 2022).

In the Villafranca project we are combining astrometry and photometry from Gaia (Prusti et al. 2016) and spectroscopy from the Galactic O-Star Spectroscopic Survey (GOSSS, Maíz Apellániz et al. 2011) and the Library of Libraries of Massive-star high-Resolution spectra (LiLiMaRlin, Maíz Apellániz et al. 2019a) to study the properties of Galactic stellar groups with OB stars. In Villafranca I (Maíz Apellániz et al. 2020a) we analyzed 16 OB groups using Gaia Data Release 2 (DR2) data and in Villafranca II (Maíz Apellániz et al. 2022) we added another ten and reanalyzed the 26 objects with Gaia Early Data Release 3 (EDR3) data. In Villafranca I we discussed the case of Villafranca O-014, the OB group associated with the North America and Pelican nebulae. Gaia DR2 data confirmed that the Bajamar star (named after its “geographical” location with respect to the North America nebula at the Bahamas), an early-type O SB2 system previously suspected to be the main ionizing source of the nebulae (Comerón & Pasquali 2005; Maíz Apellániz et al. 2016), is indeed at a similar distance as the nebula but at a location where no cluster is present (see Kuhn & Hillenbrand 2020 and Villafranca II for confirmation of the distance with Gaia EDR3). Kuhn et al. (2020) showed that the Bajamar star is a walkaway stellar system recently ejected from Villafranca O-014 NW, the most significant cluster in the region and located around Bermuda (following the geographical nomenclature). The other previously known O-type system in the vicinity, HD 199 579 (which can be dubbed the Toronto star following the geographical nomenclature), was apparently beyond the nebula using Gaia DR2 data (Kuhn et al. 2020 and Villafranca I) but Gaia EDR3 brings it within its reach (Villafranca II) when parallaxes are recalibrated according to Maíz Apellániz (2022). The Toronto star has a proper motion consistent with also being ejected from the Bermuda cluster as a walkaway system (Kuhn et al. 2020).

The North America and Pelican nebulae constitute a complex star-forming region whose pareidolic common names arise from the molecular clouds in the foreground of the H II region (Bally & Scoville 1980; Zhang et al. 2014; Kong et al. 2021). As discussed below, massive-star formation goes back at least 2 Ma but the region is still currently forming new stars, as evidenced by the study of outflows and submillimeter cores of Bally et al. (2014).

In this paper we present a scenario in which the Bajamar and Toronto stars (among other high- and intermediate-mass stars) were ejected soon after the formation of the Bermuda cluster.

2. Data and methods

In previous papers (Maíz Apellániz 2019, Villafranca I) we developed a method that uses astrometric and photometric information from Gaia to determine the distances and proper motions of OB groups. In Villafranca II we applied the method to the Bermuda cluster using Gaia EDR3 data and derived a distance d of 798 ± 6 pc, a cluster center at α = 313.10° and δ = +44.40° (J2000), and a group proper motion μ_α * ,g of −1.370 ± 0.020 mas/a and μ_δ, g of −3.054 ± 0.020 mas/a. Those values use the Gaia EDR3 astrometric calibration of Maíz Apellániz (2022) (see also Lindegren et al. 2021a,b; Maíz Apellániz et al. 2021b) and they were calculated with the correction for proper motions of Cantat-Gaudin & Brandt (2021). The calibration and the correction are also applied to the rest of the data in this paper.

In this paper we search for stellar systems that could have been ejected from the Bermuda cluster. For that purpose, we downloaded all the Gaia EDR3 sources within 15 degrees of the cluster and with G < 14 mag. For each object we applied the calibration described in the previous paragraph and computed the individual distances (see Appendix A) and relative proper motions with respect to the cluster:

$$\begin{array}{c}\hfill {\mu }_{\alpha \ast }^{\prime }={\mu }_{\alpha \ast }-{\mu }_{\alpha \ast ,\mathrm{g}},{\mu }_{\delta }^{\prime }={\mu }_{\delta }-{\mu }_{\delta ,\mathrm{g}}.\end{array}$$(1)

Based on those, we did an initial filtering of the Gaia EDR3 sources whose relative proper motions point away from the Bermuda cluster (within a tolerance angle) and whose corresponding (approximate) flight times t_f (as initially estimated from their separations with the cluster and relative proper motions) are less than 2.5 Ma. We already knew that this was the case for the Bajamar and Toronto stars (see above) and to those we added another three candidate early-B stars that passed the initial filtering (HD 195 965, HD 201 795, and HD 200 776, Tables 1 and 2) and that were already included in the ALS catalog (Pantaleoni González et al. 2021). We then traced back the trajectories of those five systems to find out if any two of them had been at the same location at the same time within the last 2.5 Ma. Doing that properly requires tracing the 3D trajectories of each pair of systems including [a] the effect of the Galactic gravitational potential and [b] the possible range of variations allowed by the uncertainties in proper motions (relatively small), distances (larger), and radial velocities v_r (a priori unknown). To trace back the trajectories in the Galaxy, we assumed a potential that can be separated in x + y (the Galactic mid-plane coordinates with x pointing toward the Galactic center and y toward the rotation direction) and z (the vertical Galactic direction toward the North). For the x + y directions we used the potential associated with a circular velocity with an exponential dependence on $r=\sqrt{x^2+{\mathit{y}}^2}$ and a characteristic length scale of 1 kpc, adjusting the values to be consistent with the distance to the Galactic center and rotation speed of the local standard of rest (LSR) of Abuter et al. (2019) and the epicycle frequency associated with the Oort constants of Li et al. (2019). For the vertical direction we used a potential derived from the mass density of Moni Bidin et al. (2012). The velocity of the Sun with respect to the LSR was taken from Schönrich et al. (2010).

As described below in each case, those initial five systems could be grouped into three different events with a relatively well determined t_f within that period frame, two each in two events and a third one with no counterpart(s). As a second step in our search for systems possibly ejected from the Bermuda cluster, we returned to the full filtered sample to detect if those ejection events had included any additional Gaia EDR3 sources. We selected the systems with distances similar to that of the Bermuda cluster and the appropriate t_f estimates and added each candidate one by one as an additional system in each ejection event to a genetic algorithm that searches through all trajectories compatible with the measured proper motions and distances (within two sigmas of each uncertainty) to find which objects could have been within 0.1 pc at the time of the event. A genetic algorithm was used at this stage because of the high-dimensionality of the problem and the need to perform tests for a large number of systems. This process yielded additional candidates that were further tested in more detail with an implementation of the TNMIN routine of Markwardt (2009) that uses the truncated Newton method to find the minimum of a function. In a first pass we did a detailed calculation that required at least one encounter with a minimum simultaneous 3D separation d_min between all systems involved in the event of less than 6 AU¹. That pass gave us our final sample of 5, 2, and 2 systems involved in each of the three ejection events, respectively. In a second pass we set up a more lenient minimum d_min of 200 AU to speed up the convergence of the search (once again, neglecting gravitational focusing) and ran Monte Carlo simulations changing the input parameters (with a two-sigma tolerance, as above) until we found 100 circumstances of each of the three events. This allowed us to calculate the statistical properties (mean, standard deviation, and correlations) of each of the input parameters and of the coordinates (α_min, δ_min) and flight times (t_min) of each event, either without weighting or using the weights from the input parameters assuming a Gaussian distribution for each one of them (Table 3).

To gain additional information about each of the nine systems involved in the three ejection events plus eight current members of the Bermuda cluster, we used the spectral types determined from GOSSS (Maíz Apellániz et al. 2011) and we derived their stellar parameters from their photometry using CHORIZOS (Maíz Apellániz 2004). Details are given in Appendix A. We also analyzed Wide-field Infrared Survey Explorer (WISE, Cutri et al. 2013) images to search for bow shocks associated with the ejected stars (Appendix B) and light curves from the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) for the nine systems to study their variability (Appendix D).

3. Results

The procedure described above yielded three ejection episodes, which we dub the Bajamar, Toronto, and HD 201 795 events, respectively, after the earliest system involved in each one of them (Table 3). They are described one by one in this section. Prior to that, we briefly analyze the motion of the (center of the) Bermuda cluster to clarify the different coordinate systems involved. We end the section with an analysis of the internal cluster motions.

3.1. The motion of the cluster center

Before we analyze the motion of the individual systems, we discuss what we know about the motion of the Bermuda cluster itself. The values above for its distance and proper motion are quite precise, as the first is known to better than 1% and the second has uncertainties of just 0.020 mas/a in each coordinate. However, its v_r is a priori unknown. It could be a significantly large value, thus implying that its distance has changed substantially in the last few Ma. To assess that possibility, we calculate the relevant velocity vector that corresponds to the measured proper motions. First, we subtract the effect in the observed cluster proper motion of the velocity of the Sun with respect to the LSR at our current location. Then, we calculate the cluster velocity in the plane of the sky with respect to the Sun LSR along Galactic longitude and latitude. Finally, we subtract the effect of the difference between the LSR at the cluster location and at the Sun’s location to obtain the velocity of the cluster with respect to its LSR in the plane of the sky. It is moving at 1 km s⁻¹ toward the Galactic West (decreasing l) and at 4 km s⁻¹ toward the Galactic North (increasing b), that is, its peculiar motion in those two directions is small, amounting to ∼8 pc in 2 Ma. Therefore, its peculiar motion in the radial direction is also likely to be of that order or lower. As 8 pc is the uncertainty in the cluster distance itself (and distance uncertainties for individual systems are larger), we can assume that in the past the cluster (and any events that took place in it) was at a distance similar to the current one. From a point of view of trajectory modeling, we will initially estimate that the ejection events took place at the same distance as the current value for the cluster but leave the value itself as a free parameter for the algorithm to determine.

A knowledge of the past trajectory of the cluster is also useful for two reasons: to calculate the ejection velocities with respect to the cluster, as that is the relevant frame to evaluate whether the ejected systems are walkaways or runaways and to estimate the conservation of linear momentum (Table 2), and to calculate the position of each event with respect to the cluster center to assess whether all originate from the same exact site or not (Table 3).

3.2. The Bajamar event

The Bajamar event involved at least five stellar systems and took place 1.611 ± 0.011 Ma (with Gaussian weights) or 1.611 ± 0.012 Ma ago (without Gaussian weights), with a total possible flight-time range in the 100 Monte Carlo simulations of 1.586–1.644 Ma. The five stellar systems we detect with a common origin are the Bajamar star, Tyc 3179-00756-1, HD 195 695, Tyc 3575-01514-1, and Tyc 3157-00918-1, with a total of at least six stars (including four massive stars) as the first one is a spectroscopic binary (Appendix C). This event is the most significant one of the three detected due to its high number of stars involved, the large total amount of mass lost by the cluster in the episode, and the inclusion of Bajamar A (the earliest and most massive star in the region, Appendix A). In the Sun’s LSR the event took place close to the current position of Tyc 3179-00756-1 (which has a very short solid cyan trajectory in Fig. 1) but if we subtract the cluster motion we find that the event took place close to its center (non-filled cyan star in the bottom panels of Fig. 1). The coordinates of the event are known with high precision (the semimajor axis of the two-sigma uncertainty ellipses in the Fig. 1 is just 2′ long) because they are the result of the intersection of five trajectories.

Fig. 1. Top left: chart in equatorial coordinates of the three ejection events in this paper, color-coded in blue (Bajamar), red (Toronto), and dark green (HD 201 795). Colored solid lines show a representative trajectory for each system in the Sun’s LSR and colored short-dashed lines the equivalent after subtracting the motion of the Bermuda cluster. The black long-dashed line shows the Galactic equator and the square indicates the location of the bottom panels. The Cygnus asterism is also plotted for reference. Top right: DSS2 image of the same field as the top left panel. Trajectories are color-coded here in cyan (Bajamar), yellow (Toronto), and green (HD 201 795). Bottom left: equivalent chart of the square region in the top panels. Stars mark the location of each event in the Sun’s LSR (filled) and after subtracting the motion of the Bermuda cluster (non-filled). The ellipses mark the corresponding two-sigma uncertainty regions of each event using Gaussian weights. The black long-dashed line shows the Galactic equator and the black crosses mark the cluster members from Villafranca II. Bottom right: GALANTE (Maíz Apellániz et al. 2021a) three-color mosaic (red: F861M+F660N, green: F665N+F515N, blue: F450N+F420N) image of the same field as the bottom-left panel. The slight curvature of the trajectories is caused by the Galactic potential.

Most of the plane-of-the-sky momentum of the event in the frame of reference of the cluster was carried away by the Bajamar star and by HD 195 965 but with a great imbalance in velocities produced by the differences in mass. Hence, the Bajamar star is still in the region (indeed, providing most of the ionizing photons responsible for the North America and Pelican nebula) while HD 195 965 is ∼6° away from the location of the event. Each of those two systems was ejected toward opposite quadrants in the reference frame of the cluster (as expected if their total linear momentum were close to zero) but when considering the four systems, there is an excess of ∼200 M_⊙km s⁻¹ toward the E. Such an excess can be corrected by tweaking the stellar masses. One possible way of doing is to reduce the mass of Bajamar A,B by 15-20% (to 85–90 M_⊙) and to increase the mass of Tyc 3157-00918-1 by 50–80% (to 5.4–6.5 M_⊙) with respect to the values derived in Appendix A. The first change is within the expected uncertainties of isochrones (especially considering variations in age) and yields a lower mass for Bajamar A that is more consistent with its spectral type. The second change would require Tyc 3157-00918-1 to be a similar-mass binary or a pre-main-sequence (PMS) star, which is a prediction that can be tested with high-resolution spectroscopy.

In this paper we only have precise tangential velocities v_t, as 3D velocities include a radial component that our modeling can only poorly constraint, as we do not have measured v_r and our Gaia EDR3 parallaxes provide only limited information about the distance differences. Considering only v_t, HD 195 965 and Tyc 3157-00918-1 are above the 30 km s⁻¹ threshold and are therefore runaway systems. The other three systems are well below the threshold and would be only walkaways unless their v_r are significantly larger than their v_t. The Gaia EDR3 distances do not appear to allow that with the only possible exception of the Bajamar star itself (which has a relatively large distance uncertainty) but a large v_r in that case would be problematic in terms of balancing momentum in the radial direction.

When we started analyzing these data, our first guess was that the Bajamar and Toronto stars were simultaneously ejected but all of our attempts with Gaia EDR3 proper motions were unable to place the two systems at the same place at the same time: at their closest approach Toronto was always SW of Bajamar. Nevertheless, we point out that the Toronto star is the brightest Gaia EDR3 source in the sample and that, as such, the multiplicative constant applied to its astrometric uncertainties is large (Maíz Apellániz 2022). If it were even larger, it would be possible to place it at the location of the event. Indeed, if we use the TGAS proper motions from Gaia DR1 (Michalik et al. 2015), a common event becomes a possibility but we will likely not have a final answer until future Gaia data releases reduce the uncertainties for systems as bright as Toronto.

3.3. The Toronto event

The Toronto event involved (at least) two stellar systems and took place 1.496 ± 0.044 Ma (with Gaussian weights) or 1.503 ± 0.051 Ma ago (without Gaussian weights), with a total possible flight-time range in the 100 Monte Carlo simulations of 1.402–1.643 Ma. Therefore, it likely happened after the Bajamar event but with Gaia EDR3 we cannot exclude that both events took place simultaneously ∼1.6 Ma ago. The two systems involved are the Toronto star and HDE 227 090, which result in three stars (Toronto is a spectroscopic binary, Williams et al. 2001) of which one or two are massive stars (Toronto B is a borderline case). In the frame of reference of the cluster, the event took place within its boundaries but the uncertainty ellipse in the lower left panel of Fig. 1 is considerably larger than that of the Bajamar event, as here there are only two trajectories involved. Most likely, the Toronto event took place to the SW of the Bajamar event but we cannot exclude that they happened at the same location. However, if that was the case, then the two events could not have taken place simultaneously given the correlations between t_min, α_min, and δ_min.

The two systems involved in the Toronto event were ejected in directions nearly perpendicular to the main direction of the previous event as defined by the Bajamar star and HD 195 965. The Toronto event is even more asymmetric, as HDE 227 090 is an intermediate-mass star ejected in the direction nearly opposite to that of Toronto (Fig. 1). Such an alignment allows for an excellent momentum conservation in the plane of sky, as evidenced in Table 2: in both right ascension and declination the total momentum including the contribution from Toronto B is within one sigma of zero. The Gaia EDR3 distance uncertainties are large (Toronto is a very bright system and HDE 227 090 has a bad renormalized unit weight error or RUWE) but the more massive system appears to have been ejected outwards with respect to the Sun and the less massive system inward, so (within the current large error bars) momentum can be balanced also in the radial direction. A corollary of this result is that the mass ratio between Toronto A,B and HDE 227 090 of 9.1 derived in Appendix A from the evolutionary masses and the mass ratio for Toronto A,B of Williams et al. (2001) agrees very well with the trajectories of the two systems assuming the absence of an additional body.

Considering its value of v_t, HDE 227 090 is a clear runaway. On the other hand, Toronto is just a walkaway in the plane of the sky. Its large uncertainty in distance does not allow us to predict its v_r but the balance of momentum in the radial direction cannot make it too large. Indeed, Williams et al. (2001) give a heliocentric center of mass radial velocity of −6.4 ± 0.4 km s⁻¹ for Toronto and in the Sun’s LSR that corresponds to +5.7 km s⁻¹. Therefore, its velocity with respect to its LSR is in the vicinity of 10 km s⁻¹, making the Toronto star a walkaway.

3.4. The HD 201 795 event

The third event took place 1.905 ± 0.037 Ma ago (with Gaussian weights) or 1.923 ± 0.062 Ma ago (without Gaussian weights), with a total possible flight-time range in the 100 Monte Carlo simulations of 1.787–2.060 Ma. We have been able to identify just two systems that participated in it, HD 201 795 and HD 200 776, both of them massive. As it happened with the previous event, uncertainties in t_min are larger because there are only two systems involved. The information on HD 200 776 is relatively limited but Dervişoǧlu et al. (2011) indicate it is an eclipsing binary, something we confirm in Appendix D, raising the number of ejected stars to three. In the frame of reference of the cluster, the HD 201 795 event happened close to the location of the previous two but the uncertainty ellipses are too large to determine an accurate separation.

One important difference with the two previous events is that linear momentum cannot be balanced in the plane of the sky, as the two systems were ejected in nearly perpendicular directions. If both systems were indeed simultaneously ejected from a 3+ body interaction, one or more stars must be missing. We searched for a third object but could not find any. One possibility is that such a third object could have already exploded as an SN. For that to happen, the interaction that led to the ejection should have happened when the third object was at least 1–2 Ma old, as otherwise there would have been no time for its evolution to lead to an SN explosion. One object that fits the separation, position angle, and approximate mass needed to balance the momentum of the event is Deneb (α Cyg), which is too bright to be included in Gaia EDR3. However, its Hipparcos distance (Maíz Apellániz et al. 2008) puts it too close to the Sun and, more importantly, its relative proper motion (van Leeuwen 2007) does not point away from the Bermuda cluster.

The value of v_t for HD 201 795 is already above the 30 km s⁻¹ threshold to be considered a runway system, independently of its v_r. HD 200 776 is below the threshold and its v_r is small (Dervişoǧlu et al. 2011) and insufficient to push the system over it, so it appears to be a walkway system. The Gaia EDR3 distance to HD 200 776 puts it beyond the North America nebula but with a large uncertainty, so there is no inconsistency with the radial velocity.

3.5. An expanding cluster

Kuhn et al. (2020) used Gaia DR2 data to detect that the Bermuda cluster (their group D) is expanding with no apparent rotation, as evidenced by their plots of tangential velocity components as a function of each coordinate (Fig. 17 in that paper) and of proper motion vectors as a function of position in the sky (Fig. 16 in that paper). They also noticed that the Bajamar and Toronto stars followed the same velocity-position trends, so “whatever phenomenon accelerated these stars to their current velocities is also responsible for the expansion of the group as a whole”, and provided flight times of ∼1.5 Ma for Bajamar and ∼1.8 Ma for Toronto. Those values are slightly different from ours, which are obtained with the more precise (thanks to the longer time baseline of EDR3 with respect to DR2) and accurate (thanks to the astrometric calibrations of Maíz Apellániz 2022) data available. In this subsection we use the new Gaia EDR3 data to reanalyze the Kuhn et al. (2020) results and in the analysis section below we discuss them in the context of how cluster orphanization takes place.

Figure 2 shoes the Gaia EDR3 relative proper motions as a function of position in the sky and Fig. 3 the behavior of ${\mu }_{\alpha \ast }^{\prime }$ and ${\mu }_{\delta }^{\prime }$ as a function of their coordinates. Those figures confirm the finding of Kuhn et al. (2020) about the expansion of the Bermuda cluster and add new information. In first place, the error bars include the effect of the covariance (significantly lower than for DR2, Maíz Apellániz et al. 2021b) and the new astrometric calibration of Maíz Apellániz (2022) to demonstrate that the expansion is real and not a measurement effect.

Fig. 2. Gaia EDR3 chart of the 144 members of the Bermuda cluster selected in Villafranca II. Symbol sizes encode G′ (corrected Gaia EDR3 magnitude), symbol colors encode $G_{\text{BP}}^{\prime }-G_{\text{RP}}^{\prime }$ (corrected Gaia EDR3 color), and arrows encode relative proper motions. Red arrows are used for the proper motions of the 17 brightest systems in G′. A cross indicates the used cluster center.

In second place, there is a significant behavior as a function of G′ magnitude that, though affected by extinction, is highly correlated with mass (see Appendix A and Fig. A.1). More specifically, we divided the 144 Bermuda members selected in Villafranca II into two groups by magnitude: the 17 brightest systems in one and the 127 faint ones in the other. When we do that, the second group is clearly more correlated in the left plot of Fig. 3 (expansion in right ascension) and the effect is even more pronounced in the right plot (expansion in declination), where the second group is still significantly correlated while the bright systems show essentially no correlation (r = 0.01). In the declination plot the effect may be described as the presence of two faint populations, a dominant one that follows the ${\mu }_{\delta }^{\prime }$ versus Δδ trend and a secondary one in the upper left quadrant. When looking at the bright systems, no trend is visible at all and the 17 points show an uncorrelated scatter about the center of the plot. A similar effect is seen in Fig. 2, where the majority of the black arrows point away from the center and their size increase with their distance to it while red arrows do not show a clear pattern and instead appear to point in random directions. We have applied a weighted linear regression with an outlier-detection algorithm to fit the data for the faint systems in both panels of Fig. 3 to obtain the expansion timescale in each case. In right ascension we obtain 1.852 ± 0.077 Ma (104 systems) and in declination we obtain 1.627 ± 0.095 Ma (107 systems). Those values are discussed in the analysis section below.

Fig. 3. Relative proper motions in right ascension (left) and declination (right) for the 144 members of the Bermuda cluster selected in Villafranca II as a function of the difference in each coordinate with respect to the cluster center. Red symbols are used for the 17 brightest systems in G′ and black symbols for the 127 fainter systems. The corresponding separation and velocity scales are given in the top and right axes, respectively. The correlation factors for the bright and faint systems are also given. The region in gray shows the possible linear fits to the faint systems in each case after rejecting outliers to leave 104 objects (left) and 107 (right).

In third place, the Bermuda cluster appears to have been hollowed out in Fig. A.1. There is only one system close to its center and is a bright one, i.e a member of a group that as we have seen appears to follow a different pattern with respect to the (majority of) faint systems. This pattern is the opposite of what would be expected for a centrally concentrated cluster or even for one with a spherically symmetric density distribution with a flat value near its center. Therefore, it is another piece of evidence in favor of the expansion of the cluster.

Kuhn et al. (2020) proposed three possible mechanisms for the expansion of the Bermuda cluster: tidal fields, molecular gas dispersion, and initial velocity gradient in their natal cloud. In the analysis section below we propose a fourth possibility, namely the mass loss associated with the three ejection events in this paper.

4. Analysis

The discovery of the multiple stellar ejections from the Bajamar cluster leads to several questions, which are analyzed here. We start with the confidence on our results, we follow with the effect that the ejections have on the mass function, and we then discuss how cluster orphanization happens. Finally, we discuss other issues related to the discovery and with possible ways to extend this work.

4.1. On the reality of the ejection events

In principle, it could be possible that what we are seeing is just a coincidence and that the ejection events did not really occur, with the trajectories just happening to overlap. To explore such a possibility we start with the Bajamar event. The first obvious point that needs to be made is that it is not simply that the trajectories overlap but that they coincide at the same position of the sky at the same time. This is relevant because the majority of the systems at the range of distances we are sampling have small relative velocities with respect to their LSR that are not sufficient to cross several tens of pc in less than 2 Ma, as for example, HD 195 965 and Tyc 3157-00918-1 have done. With that in mind, the probability of back-tracing five trajectories within the sample volume (a truncated cone with a radius of ∼200 pc and a depth of ∼100 pc) to the same location within ∼200 AU or less is extremely low. Furthermore, all five systems are of early-type (O, B, or early-A), including the earliest O-type object in the region (or indeed within 1 kpc of the Sun altogether), another possible O star, and an early-B star. At a distance of ∼800 pc, the sample is dominated by F and early-G stars with G fainter than 12 but we see none of those participating. Indeed, despite sampling down to G = 14, the faintest star detected in the Bajamar event, Tyc 3575-01514-1, has G = 11.0, so the ejected systems are not randomly drawn from the sample at all. Additionally, the fact that only minor adjustments to the evolutionary masses are needed to balance momentum is hard to explain without the alleged simultaneous multiple ejections being a real physical event. Finally, despite allowing for a relative large search area in the sky for the encounter, the location of the event is a region of less than 1 pc in size very close to the cluster center (once its motion has been subtracted). Therefore, while it is true that our algorithm is searching essentially in 2D space (as we do not have much information about depth), the above circumstances make it extremely unlikely that the Gaia EDR3 data allowing all five systems within such a small volume at the same point in time is just a coincidence.

For the Toronto event we can make some of the same arguments as for the Bajamar event but in this case only two systems are involved and their astrometry has larger uncertainties. As a result, the location of the event, though still very close to the cluster center, has a larger uncertainty. Nevertheless, in this case the momentum balance is even more impressive than in the previous one. What we have found are two systems ejected in opposite directions (within the uncertainties) and with very different masses (q = 9.1) but with velocities that compensate for the mass ratio perfectly (within one sigma) to conserve momentum, a feat difficult to accomplish with just two systems. Furthermore, the age of the event is very similar to that of the Bajamar event, a connection that points toward a large fraction of the cluster forming simultaneously if ejections take place preferentially within a short time after its birth. Therefore, the case for the reality of the Toronto event is still strong. The one point that could change in the future with better astrometry would be the relationship between the Bajamar and Toronto events, that is, their degree of spatial and temporal simultaneity.

The HD 201 795 event is the one for which the certainty is lower. As for the Toronto event, there are only two systems and that leads to a larger uncertainty on its location (though the data are still compatible with the event taking place near the cluster center). More importantly, momentum is clearly imbalanced so either a third body is missing or the event did not really happen. However, we should also point out that this is the only event for which all the detected systems (two in this case) are bright and previously well identified as being of early type, for example both were present in the ALS catalog and both are seven-magnitude objects in G. As a reference, they are among the 400 brightest (in G) systems within 15 degrees of the Bermuda cluster with measured Gaia EDR3 parallaxes between 0.8 mas and 2.0 mas. The importance of the HD 201 795 event is that, if real, it extends the time span of the ejection events ∼300 ka in the past from the Bajamar event which, if the ejections really take place a short time after the massive stars are formed, implies that there are several star-forming events separated by several hundreds of ka. We return to that point below but here we finish by stating that the confirmation of the reality of the HD 201 795 event is the most pressing issue that still needs to be resolved regarding the orphanization of the Bermuda cluster.

4.2. The IMF and the PDMF

The preferential ejection of massive stars necessarily modifies the mass function of the cluster. We cannot carry a detailed analysis of the mass function because the presence of strong differential extinction in the Bermuda cluster (Appendix A) would require the use of IR information to characterize it and to estimate the completeness of the sample for low-mass stars. For that reason, we will center our analysis on the intermediate- and high-mass population.

We use the data in Appendix A to determine the number of stars in the 12–100 M_⊙ and 4–12 M_⊙ ranges before and after the ejection events. Of the nine ejected stars listed in Table A.2, six are in the first range and just one on the second. Adding Bajamar B and Toronto B changes those numbers to seven and two, respectively. Of the current eight cluster members, two are in the 12–100 M_⊙ range and six are in the 4–12 M_⊙. The number of stars in the first range is very likely to be complete. As for the second range, converting Tyc 3157-00918-1 into a binary does not change the accounting and the eclipsing companion to HD 200 776 is unlikely to be as massive as 4 M_⊙. Regarding unaccounted current cluster members, a look at Fig. A.1 indicates the existence of several stars in the 2–4 M_⊙ range and it is possible that one or two of them may actually be over 4 M_⊙. In summary, before the ejections the ratio between the two ranges was 9:8 (or maybe 9:9 or 9:10 if there are missing stars) and the current value after the ejections is 2:6 (possibly 2:7 or 2:8). The total ejected mass from Table A.2 plus Bajamar B, Toronto B, and HD 200 776 is ∼220 M_⊙. Considering the possible overestimation of the Bajamar mass discussed above, the real value may be slightly smaller but we can establish a lower limit of 200 M_⊙.

We compared this with a Kroupa IMF, which is a multi-part power law with slopes of −0.3 (M < 0.08 M_⊙), −1.3 (0.08 M_⊙ < M < 0.50 M_⊙), and −2.3 (0.50 M_⊙ < M) (Kroupa 2001). Under that assumption, if there are eight stars in the 4–12 M_⊙, one expects 2.4 stars in the 12–100 M_⊙ range (changing the number to ten stars in the first range raises the second number to 3.0). The observed number of 9 (before ejections) is over three sigmas away from the expected value. Therefore, we conclude that the IMF of the Bermuda cluster was significantly top-heavy, with a measured value of the intermediate/massive-star slope of ∼ − 1.3 instead of the Kroupa value of −2.3. On the other hand, in the current configuration of the cluster there are six stars in the 4–12 M_⊙ range and in that case a Kroupa mass function predicts 1.8 stars with masses of 12–100 M_⊙, which is in good agreement with the observed value of two. This leads to the interesting conclusion that cluster orphanization can transform a top-heavy IMF into a Kroupa-like PDMF. This discovery is even more relevant when we consider that the Bermuda cluster is the only known case of a stellar group where stars as massive as Bajamar A have been recently formed within 1 kpc of the Sun and, therefore, is one of the clusters that we can potentially study better: what we see here may be hiding somewhere else.

Another consequence of the assumption of a Kroupa IMF is that if there are eight stars in the 4–12 M_⊙, one expects ∼655 cluster members in the 0.08–100 M_⊙ range. Kuhn et al. (2020) identified 235 members of Bermuda cluster but not reaching down to brown dwarfs so the numbers are consistent within a factor of two. Such a cluster would have a total initial mass of M_clus∼376 M_⊙. The real Bermuda cluster may have had a slightly different initial mass, higher by ∼150 M_⊙ for being top-heavy and possibly lower if the Kuhn et al. (2020) measurement is more complete than expected and the cluster was also initially top-heavy when comparing intermediate-mass stars and low-mass ones. A conservative assumption is that the Bermuda cluster formed 400–500 M_⊙ in stars, making it a moderately low-mass cluster. If such a cluster loses ∼200 M_⊙ soon after its formation, when an additional mass loss in the form of gas is also likely to happen, it could be sufficient to unbind it. We explore that possibility in the next subsection.

In Villafranca II we initially proposed that the Bajamar star had formed in near isolation. After its characterization as a walkaway star by Kuhn et al. (2020) we clarified that by near isolation we included the possibility of the star being born in a low-mass cluster, something that was considered highly unlikely in the Weidner & Kroupa (2006) scenario but that the results in this paper confirm. The hidden variable turned out to be the multiple ejections of (preferentially) massive stars significantly changing the PDMF with respect to the IMF.

4.3. How cluster orphanization happens

The results presented in this paper make the Bermuda cluster a Rosetta stone for understanding massive-star formation and cluster evolution for two reasons: there is no other stellar cluster that has recently formed stars as massive as Bajamar A within 1 kpc of the Sun and it is extremely young, with episodes younger than 2 Ma. In the previous subsection we analyzed the consequences that the ejection of the most massive stars of a young cluster has on the PDMF; here we discuss what the process tells us about the formation of massive stars and the evolution of stellar clusters in their initial stages.

Numerical simulations tell us that the critical parameter to produce dynamical ejections is the cluster radius: ejections are more common in compact clusters, as expected from the increased likelihood of close encounters at high densities (Oh & Kroupa 2016). Also, for compact clusters with masses lower than 1000 M_⊙ simulations indicate that it is possible to eject all stars more massive than 17.5 M_⊙ and produce an orphan cluster (Oh et al. 2015). Those papers assume a Kroupa IMF but, as we have seen above, the Bermuda cluster has a top-heavy IMF and that increases the possibility of ejections given the existence a higher fraction of massive stars. Therefore, the events described in this paper are likely to have happened in an environment quite different to that of the present Bermuda cluster, when the cluster was highly compact. In this respect, above we derived expansion timescales for the Bermuda cluster in the 1.6–1.9 Ma range, values that are compatible with the age of the ejection events². As we have seen above, the ejected mass may be sufficient to produce such an expansion. All of this points toward a common origin for the ejection events and the cluster expansion at a time when the Bermuda cluster was highly compact.

Nevertheless, the above scenario does not explain two issues: the significant age discrepancy between the HD 201 795 event and the other two and the fact that the most massive stars still in the cluster (along with a few low-mass ones) do not partake on the cluster expansion. One possible explanation is that clusters do not form monolithically but instead do so in a “conveyor-belt” fashion (Longmore et al. 2014; Krumholz & McKee 2020), with two molecular filaments flowing in opposite directions and colliding near-continuously during timescales (> 1 Ma) longer than those required to form massive stars. This allows large numbers of stars to be formed but the timescale poses a problem: once an O star is born, it will ionize its surroundings and hamper the subsequent formation of O stars. In the context of the Orion Nebula Cluster, Kroupa et al. (2018) suggested a way around: if O stars are preferentially ejected as runaways soon after their births, new generations of star formation can take place, which brings us back to the Bermuda cluster. We have seen that it has been a highly efficient location for massive star ejections, likely because of a very small radius, and this could have been the cause for at least two generations of massive-star formation: the first one that produced HD 201 795 and HD 200 776 and a second one that produced Bajamar, Toronto, and most of the rest of the stars. Furthermore, after ejecting all the stars form the three events, star formation could have kept going and produce some additional stars. As those would not have been implicated in the events, they would not follow the cluster expansion pattern and instead form the contaminant population in Figs. 2 and 3. The separation of the star formation in the cluster into different episodes separated by hundreds of ka also simplifies the energy balance of the cluster expansion: the major episode would have been the one that involved Bajamar and Toronto but the more recent one would have added stellar mass to the cluster a posteriori, making the Bermuda cluster even less massive at the time of the ejections. The high rate of stellar ejections would also explain the anomalous top-heavy IMF: as O stars have not been left around to ionize the infalling molecular gas, massive-star formation has continued for a longer period of time than if they had remained in the cluster. In summary, the available data are consistent with star formation taking place in the Bermuda cluster in a conveyor-belt fashion with massive stars being ejected and in that way suppressing their negative feedback and leading to a top-heavy IMF.

As we discussed in Villafranca I, there are two views on the formation of OB associations: the traditional (or nurture) view of Lada & Lada (2003) that states that all stars are born in clusters but the feedback from protostellar outflows, UV radiation, and stellar winds sweeps the leftover gas and unbinds the stars by lowering the absolute value of their potential energy; and the alternative (or nature) view of Elmegreen (2010), Ward et al. (2020) that states that star formation is a hierarchical process that happens at very different length scales, leading to stars being born in both bound (clusters) and unbound (associations) stellar systems. The analysis of Krumholz & McKee (2020) suggests that clusters are formed under the special circumstances of the conveyor-belt model when the configuration is such that the molecular gas in two opposing filaments keeps feeding the same location of star formation for a significant amount of time. In most circumstances this does not happen and we end up with multiple sites of star formation in a large cloud. Also, the full mass of a single cloud is never assembled at a single point in time. Recent analyses of the expansion of OB associations have produced mixed results (Wright 2020) but, in general, it is not true that all show a consistent isotropic expansion, tipping the balance toward the hierarchical model of star formation. The next question is how the Bermuda cluster fits in this context. Interestingly, this is a system that appears to have been bound at one point and where the main possible source of energy injection into the leftover gas (O stars) was lost early in the process, thus facilitating the continuation of the star formation process through the conveyor-belt mechanism. However, the cluster eventually lost a significant amount of mass in the form of stars (and, at some point, also gas) and has ended up unbound, so its future is an expanding OB association. We would like to know how common is such a process. We know that the Bermuda cluster is not a run-of-the-mill system, as stars as massive as Bajamar A are seldom found in star-forming regions but in this respect, as in most aspects related to the top of the IMF, the sample is scarce. Further analyses with Gaia of similar systems, looking for evidence of similar ejections, are needed.

A final point about the process of cluster orphanization is that, as already mentioned, we have been unable to detect any ejected star below 3 M_⊙ despite having data of sufficient quality to detect systems down to ∼1 M_⊙ with moderate extinction. A scenario that explains that is based on the conveyor-belt model described above. Star formation in the Bermuda cluster took place during ∼1 Ma in a series of events spaced by hundreds of ka. In each event a combination of low-, intermediate-, and (in some of them) high-mass star formation was started but before the next event happened, only the most massive systems had time to reach the zero age main sequence (ZAMS) and the rest were still on the PMS. In the time between events, the newly formed stars were in unstable dynamical configurations that eventually led to the three ejection events described in this paper. However, any low-mass stars would have been significantly larger in size (especially if still on the early phase of their Hayashi tracks) at that point and, if involved in such a dynamical interaction, they could have been disrupted by tidal forces. This would explain the top-heavy mass distribution of the ejected objects and would indicate that ejections take place very early in the life of the cluster.

4.4. Other issues

There are several early-type O stars among the previously confirmed runaway/walkaway stars. Gvaramadze & Gualandris (2011) mention BD +43 3654 (O4 If, Maíz Apellániz et al. 2016 and λ Cep (O6.5 I(n)fp, Sota et al. 2011) and Maíz Apellániz et al. (2018) mention ALS 4962 (an uncertain ON5 Ifp) and ALS 11244 (O4.5 III(n)(fc)p). However, those are later than the three stars classified as O2 If*/WN5 or O2 If*/WN6 ejected from Westerlund 2 (Villafranca O-004): THA 35-II-42 (= WR 21a), SS 215 (= WR 20aa), and WR 20c (Roman-Lopes et al. 2011; Villafranca I). Outside the Milky Way, the most massive runaway known is VFTS 682 (Renzo et al. 2019b). To that select club of massive runaways/walkaways we add the Bajamar star.

The next question is at what stage ejections happen. In the numerical simulations of Oh & Kroupa (2016) most ejections take place before an age of 1 Ma but some happen as late as 3 Ma. However, those simulations are for a monolithic cluster with 3000 M_⊙, so they do not necessarily apply here. If the conveyor-belt scenario above is correct, then most ejections may happen even earlier than 1 Ma, right after the first massive stars are formed and the cluster is in a highly dynamical state caused by the chaotic interactions at its core (Bate et al. 2003). This would be consistent with the agreement between the two CHORIZOS results for Bajamar A in Appendix A, which imply that the age of the system is not too different from that of the 1.78 Ma isochrone. On the other hand, if the reason why the HD 201 795 event does not have a balanced momentum is because the third body has already exploded as an SN, then the ejection must have taken place at least 1 Ma (and possibly more) after the formation of the progenitor. A further analysis of the system may provide answers to the question but, at the very least, our study sets a lower limit for the age of the first massive stars in the Bermuda cluster (1.9 Ma ago) and for the age of the first O stars (1.6 Ma ago).

Other stars in the region may have been ejected from the North America nebula. The Cygnus loop or Veil nebula is a well studied supernova remnant (SNR) at a distance of 735 ± 25 pc (Fesen et al. 2018) to the south of the North America nebula just outside the top panels in Fig. 1. We want to know if the progenitor may have been ejected from the North America nebula. A velocity of ∼100 km s⁻¹ would be needed if the age of the event is 2.0 Ma but a lower value would suffice if it were older (thus placing the north of the first massive stars in the region even further into the past). Fang et al. (2017) already modeled the morphology of the Cygnus loop as produced by a runaway star with characteristics similar to those described here but we must point out that the direction from the nebula to the SNR is the wrong one to balance the momentum of the HD 201 795 event, so if there is a connection then it must have originated in a separate event. One prediction we can make is that as the distance is shorter than the value for the Bermuda cluster, chances are that if the progenitor was ejected from the North America nebula, then it was originally moving toward us. In that case, the southern edge should be closer than the northern edge.

As discussed in Villafranca I and II, other clusters may have experienced orphanization. Villafranca O-012 S (Haffner 18) is probably the object that deserves to be analyzed next, as its most massive stars left the cluster ∼400 ka ago. Villafranca O-023 (Orion nebula cluster), arguably the best-studied star-forming region that has produced O stars, is another complex case in which objects were ejected as early as 2.5 Ma ago (Blaauw & Morgan 1953; Hoogerwerf et al. 2000) and as late as in historical times (Bally et al. 2020; Maíz Apellániz et al. 2021c). In that case the question should not be whether the most massive star ever formed there is still in the cluster (it apparently is, θ¹ Ori Ca) but whether at some point the most massive stars were ejected and left the cluster temporarily orphaned (apparently that happened with the event from 2.5 Ma ago). If anything is clear, it is that the formation of massive stars and their relationship to stellar clusters is so complex that each case needs to be analyzed in detail in order to understand it.

4.5. Future work

There are several avenues that can be explored to confirm and improve upon the results presented here. The most obvious is the measurement of the radial velocities, which requires multiple epochs to discard the existence of spectroscopic binarity (if single) or the measurement of the systemic velocity γ (if multiple) and a consistent methodology to avoid the offsets of several km/s that plague the values for OB stars (Trigueros Páez et al. 2021). Second on the list is determining spectral classifications for those objects without them and, possibly, detailed spectroscopic modeling to refine the stellar parameters. Third is the derivation of the spectroscopic orbits for objects such as Bajamar. Finally, the improved astrometry from future Gaia data releases should place more astringent values on the properties of the analyzed stars.

Acknowledgments

We thank the referee, John Bally, for his suggestions to improve the paper. J. M. A. thanks M. A. Kuhn for a conversation that gave him the idea to search for stars escaping from the North America nebula and M. Cerviño for useful comments. J. M. A. and M. P. G. acknowledge support from the Spanish Government Ministerio de Ciencia e Innovación through grant PGC2018-095 049-B-C22. R. H. B. acknowledges support from ANID FONDECYT Regular Project 1 211 903 and the ESAC visitors program. M. W. acknowledges support from the Spanish Government Ministerio de Ciencia e Innovación (MICI/FEDER, UE) through grant RTI2018-095 076-B-C21, and from the Instituto de Ciencias del Cosmos de la Univ. de Barcelona (ICCUB, Unidad de Excelencia “María de Maeztu”) through grant CEX2019-000 918-M. 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. The Gaia data are processed with the computer resources at Mare Nostrum and the technical support provided by BSC-CNS. This publication makes use of data products from the Wide-field Infrared Survey Explorer (WISE), which is a joint project of the University of California, Los Angeles; and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration (NASA), where the word “National” refers to the United States of America. This paper includes data collected by the Transiting Exoplanet Survey Satellite (TESS) mission. Funding for the TESS mission is provided by the NASA’s Science Mission Directorate. This research has made extensive use of the SIMBAD and VizieR databases, operated at CDS, Strasbourg, France.

References

Appendix A: CHORIZOS analysis

A.1. Data and methods

To derive the properties of the systems in this paper from their photometry we use the Bayesian code CHORIZOS (Maíz Apellániz 2004).

Photometric data. We use six optical photometric bands, three from Gaia DR2 and three from Gaia EDR3, G_BP, DR2 + G_DR2 + G_RP, DR2 + G_BP, EDR3 + G_EDR3 + G_RP, EDR3, applying the corrections and calibrations of Weiler et al. (in preparation). As discussed there, the effective passbands have changed slightly from DR2 to EDR3 and additional information can be obtained by combining the two data releases rather than using just one of them. To do that accurately, one needs to apply several corrections, divide the G_BP, DR2 and G_RP, DR2 into two G_DR2 magnitude ranges at a value of 10.87, and do a similar division for the DR2 and EDR3 G photometry at a magnitude of 13.0. We note that the new DR2 photometric calibration is similar but slightly different from that of Maíz Apellániz & Weiler (2018), as the EDR3 reanalysis has allowed for a better understanding of the DR2 issues. We also use the three 2MASS near-infrared (NIR) bands J_2M + H_2M + K_2M (Skrutskie et al. 2006) with the Maíz Apellániz & Pantaleoni González (2018) calibration. To check for the existence of mid-infrared (MIR) excesses we have also downloaded their WISE W1 + W2 + W3 + W4 photometry (Cutri et al. 2013) but we have not included it in the fit itself. As listed in Tables A.1 and A.2, for one system we have discarded the G_EDR3 photometry because it was incompatible with the other five optical bands and for another three we have eliminated one or two of the 2MASS bands because their IR excesses extended into the NIR.

Models. We use two grids based on the T_eff-luminosity class (LC) grid with solar metallicity from Maíz Apellániz (2013b). One is the full grid in which T_eff and LC are the two intrinsic independent parameters and for the other one we force the result to be in the 1.78 Ma solar-metallicity Geneva isochrone with no rotation (Lejeune & Schaerer 2001), for which the only intrinsic independent parameter is the stellar initial mass, m_i. The parameter range that is investigated is covered by the TLUSTY (hot stars, Lanz & Hubeny 2003, 2007) and Munari (cool stars, Munari et al. 2005) parts of the grid, so data from both sources are used. In particular, we note that the NIR colors from TLUSTY models for late/mid-B stars are incorrect by up to several hundredths of a magnitude, so in that region of the spectrum the Maíz Apellániz (2013b) grid uses the Munari models for all hot stars.

Distances. For the systems that are still cluster members we use the distance to the Bermuda cluster from Villafranca II (798 pc). For the walkaways/runaways we use the Gaia EDR3 individual parallaxes with the calibration of Maíz Apellániz (2022) and use the prior of Maíz Apellániz (2001, 2005) with the parameters of Maíz Apellániz et al. (2008). In particular, for the walkaways we assume they are part of the disk population and for the runaways that they are part of the halo or runaway population.

Other parameters. For the full T_eff-LC grid runs we fix the distance and fit the other four parameters (T_eff, LC, amount of extinction E(4405 − 5495), and type of extinction³R₅₄₉₅) unless the star has a known spectral type, in which case we also fix T_eff. To establish T_eff from the spectral type we use a version of the Martins et al. (2005) calibration extended to later spectral types and adapted to the subtypes changes of Sota et al. (2011, 2014), Maíz Apellániz et al. (2016). For the isochrone grid fits we also fix the distance and fit the other three parameters (m_i, E(4405 − 5495), and R₅₄₉₅). For all runs we use the family of extinction laws of Maíz Apellániz et al. (2014).

On using two grids. The reasons for using two grids (full T_eff-LC and isochrone) is that each one has its pros and cons and that two independent results allow for additional checks. The most accurate results are those for the T_eff-LC grid when T_eff is fixed from the spectral type. In that case if a system is a hidden multiple with similar spectra types the effect is to make the system of a higher LC without altering other parameters, as colors for OB stars have a weak dependence on luminosity. On the other hand, if T_eff is not fixed because no spectral type is available, the resulting uncertainties can be rather large (see e.g. Tyc 3179-00756-1 in Table A.1) as the filter combination used provides some T_eff discrimination but not much (Maíz Apellániz et al. in preparation). The situation becomes worse for stars with high extinction and for that reason cluster systems are fitted only using the isochrone grid (Table A.2) and not the full T_eff-LC grid (Table A.1). The isochrone grid fits works very well if the age is well known and if there are no vertical shifts in magnitude caused by distance errors or hidden companions. For that reason, for the isochrone results for the Bajamar and Toronto stars in Table A.2 we subtracted the effect of their companions by adding 0.18 mag and 0.05 mag, respectively, to all magnitudes, with the values estimated from the available spectroscopy and the expected differences in spectral types. Using an isochrone of the wrong age can create significant biases in some parameters and the resulting uncertainties can be underestimated, as it is the case for m_i, T_eff, and log L in Table A.2.

A.2. Stellar results

We analyze the CHORIZOS results for the systems in Tables A.1 and A.2 dividing them by the event they belong to and finishing with the ones that are still cluster members. The extinction characteristics are discussed in the next subsection.

Bajamar event. The Bajamar star, the most massive system in this paper, shows an excellent consistency in T_eff between the two runs, indicating that its age cannot be too different from the 1.78 Ma of the isochrone, and a value of 1.13 for the LC, indicating that it has already evolved from the ZAMS. The differences in G_abs between the two runs can be adscribed to the second one referring only to the A component (of spectral type O3.5 III(f*)). In conjunction with the mass ratio of 2.2±0.5 derived below, these results indicate that the Bajamar star is a binary with masses of ∼73 M_⊙ and 33±8 M_⊙, the latter within the range expected for an unevolved O8 dwarf⁴. The first run results for Tyc 3179-00756-1 have relatively large uncertainties but those from the second run point out toward the object being an O9.7/B0 dwarf if it is a single star of age younger than 2 Ma. This could be a third O star ejected from the Bermuda cluster, something that should be confirmed spectroscopically. The two runs show similar results for HD 195 965, albeit slightly hotter for the isochrone run, which translates into a higher log L due to the larger bolometric correction. The last two systems involved in the event, Tyc 3575-01514-1 and Tyc 3157-00918-1, have results compatible with being late-B or early-A dwarfs. Tyc 3575-01514-1 has a MIR excess in the W3 and W4 WISE bands. The results for the two runs for Tyc 3157-00918-1 are somewhat discrepant, with those of the first run indicating a lower T_eff. One possible way to reconcile them would be if the system is actually composed of two similar early-A stars, as that configuration will lead to the T_eff of the first run but a higher luminosity that would be confused in the second run with a higher T_eff. Alternatively, it could be a PMS star of A type. As discussed in the text, in either case the total mass of Tyc 3157-00918-1 would be higher than the value in Table A.2, which assumes a single MS star.

Toronto event. The Toronto star shows a good agreement in T_eff between the two runs, with a slightly higher value in the second one but within the uncertainties expected from the T_eff-spectral type calibration (Maíz Apellániz et al. 2014), also pointing toward an age close to that of the isochrone. Assuming the mass ratio of 4 from Williams et al. (2001) leads to a system composed of a ∼33 M_⊙ and an ∼8 M_⊙ stars, the second most massive of the ones analyzed in this paper. HDE 227 090 is an intermediate-mass star with a lower value of T_eff using the T_eff-LC grid than using the isochrone grid. Its predicted spectral type is mid/late-B. We note that the G_EDR3 magnitude of HDE 227 090 is incompatible with the other five Gaia photometric points and has been excluded from the analysis.

HD 201 795 event. The two systems involved in this event have results in the two runs that are consistent with both being early-B dwarfs. In the case of HD 201 795 the consistency between the two runs is excellent. Those of HD 200 776 are slightly discrepant, with the Table A.2 results pointing toward higher T_eff values than those expected from its spectral type. The discrepancy could be resolved in two ways: if the star has a spectral type of B0.2-B0.5 V instead of the published one (in this case we do not currently have GOSSS or LiLiMaRlin data) or if the star turns out to be closer to us by 50-100 pc, something that is within the uncertainties of the current Gaia EDR3 distance and indeed is suggested by the radial velocity measurement of Dervişoǧlu et al. (2011) assuming that the two systems were close to the current distance to the Bermuda cluster at the time of the event. The mass ratio of 5.2 derived by Dervişoǧlu et al. (2011) indicates that the companion is an intermediate-mass object of ∼2.6 M_⊙ and, therefore, has only a minor contribution to the combined light output (Appendix D and Fig. D.1).

Cluster members. The brightest cluster members of the Bermuda cluster selected in Villafranca II have been analyzed using the isochrone grid. They all have significant extinctions (see below) and have not been studied in much detail due to the recent identification of the cluster. They appear to be two stars with masses of 12-25 M_⊙, three with masses of 8-12 M_⊙, and another three with masses of 4-8 M_⊙. Three of them, 2MASS J20531694+4424298, 2MASS J20512895+4404230, and V2022 Cyg have IR excesses that reach into the NIR and for another two, 2MASS J20522503+4437562 and Tyc 3179-00733-1, the excess is visible but only in the W3 and W4 WISE bands. The most interesting object among the cluster members is 2MASS J20531694+4424298, the most extinguished target of the eight and with a photometry consistent with being a late-type O star. However, such an identification needs a spectroscopic confirmation, especially considering the anomalous spectral energy distribution that this object has in the IR. A second object, 2MASS J20522503+4437562 is another potential O9.7/B0 star that deserves spectroscopic observations to confirm such spectral type.

A.3. Extinction results

The systems analyzed in Tables A.1 and A.2 show a wide range of extinctions (see also Fig. A.1). Those still in the cluster or in its neighborhood (the Bajamar star and Tyc 3179-00756-1) have large or very large extinctions while those that have been ejected at large separations have low values of E(4405 − 5495) around 0.2 mag, which are typical for objects at that distance in this region of the Milky Way (see Fig. 8 in Maíz Apellániz & Barbá 2018).

Fig. A.1. Color-magnitude diagram with the extinction-corrected values (filled color circles) from Table A.2 and extinction-uncorrected ones (unfilled color circles). Some systems are labeled to allow for an easier correspondence with Table A.2. Different colors are used for the Bajamar (red), Toronto (blue), and HD 201 795 (dark green) events and for the eight brightest members of the Bermuda cluster (magenta). The Bajamar B and Toronto B components are not plotted but are discussed in the text. Unfilled black circles show the extinction-uncorrected measurements for the rest of the Bermuda cluster members selected in Villafranca II. The two black lines are the 1.78 Ma solar-metallicity Geneva isochrone with no rotation (a) without extinction and (b) with E(4405 − 5495) = 1.0 mag and R5495 = 3.0. Magnitudes are from Gaia EDR3.

The Bajamar and Toronto stars had been previously analyzed with different data sets in Maíz Apellániz et al. (2021d) and Maíz Apellániz & Barbá (2018), respectively. The new results are in reasonable agreement with the previous ones but uncertainties are lower now thanks to the better quality and improved calibration of the Gaia photometry. Nevertheless, we should point out that while the ${\chi }_{\text{red}}^2$ values are close to 1.0 for low- and intermediate-extinction objects, the three most extinguished systems have higher values, with the worst case being the Bajamar star itself, which is the most extinguished of all. As pointed out in Maíz Apellániz et al. (2021d), such an increase in ${\chi }_{\text{red}}^2$ with extinction is a sign that the family of extinction laws of Maíz Apellániz et al. (2014) is starting to lose validity for objects as extinguished as the Bajamar star (likely as a result of the wrong slope in the NIR, see Maíz Apellániz et al. 2020b) but we currently do not have a better alternative, as the other popular choices of Cardelli et al. (1989) and Fitzpatrick (1999) fare even worse (Maíz Apellániz & Barbá 2018).

The results in this paper also confirm other conclusions of Maíz Apellániz & Barbá (2018) and Maíz Apellániz et al. (2021d). Extinction can vary greatly within a cluster not only in amount but also in type of extinction, so using a single value of R₅₄₉₅ can lead to significant biases. In the case of the North America nebula (excluding runaway systems already far from it) we see values of R₅₄₉₅ that range from ∼2.7 to ∼3.6. Furthermore, the changes are clearly associated with the nebular morphology. Objects behind dense molecular clouds (e.g. the Bajamar star or Tyc 3179-00756-1) tend to have large values of E(4405 − 5495) and small ones of R₅₄₉₅ while objects immersed in the H II region (e.g. the Toronto star) tend to have small values of E(4405 − 5495) and large ones of R₅₄₉₅. As we proposed in Maíz Apellániz & Barbá (2018), extinction in H II regions appears to be produced mostly by large dust grains due to the selective destruction of small grains, which are instead abundant in molecular clouds.

Appendix B: Bow shocks

We have searched for bow shocks created by the walkaway and runaway systems in their WISE W3+W4 images, as such structures can be seen in the WISE images if the conditions are appropriate (Maíz Apellániz et al. 2018). The Bajamar star and Tyc 3179-00756-1 (which is relatively nearby in the plane of the sky) are immersed in a bright MIR nebulosity, likely caused by the warm dust on the far side of the Atlantic Ocean molecular cloud heated by them, that does not allow for any faint structure such as a bow shock to be detected. On the other hand, the other four walkaways/runaways that are known to be of spectral type O or early-B do show bow-shock-like structures (Fig. B.1). The O star (Toronto) displays the largest one, followed by the earliest B-type, HD 195 965, as expected from their wind strengths. The bow shocks are better aligned with the relative proper motions than with the absolute ones (especially so for the Toronto star and HD 200 776), as expected if the ISM in the North America nebula and its surroundings have similar proper motions to that of the Bermuda cluster. These detections and their orientation confirm the likelihood of the ejections of these four systems from the Bermuda cluster.

Fig. B.1. WISE W4+W3+W2 red-green-blue (RGB) mosaics for four walkaway and runaways systems. Each field is 20′×20′ (4.7 pc × 4.7 pc at a distance of 798 pc) and is oriented with north toward the top and east toward the left. In each mosaic the runaway candidate is at the center and the arrows show the absolute (green) and the relative-to-the cluster proper motions (yellow).

Appendix C: On the SB2 nature of the Bajamar star

Maíz Apellániz et al. (2016) identified the Bajamar star as an SB2 and indicated that some epochs showed velocity differences of ∼300 km/s between the He I and He II lines. No SB2 orbit has been published at this time, most likely because of the difficulty of studying this system. The spectra is dominated by the H and He II lines but the extremely high extinction makes the blue-violet region difficult to access. The motion of the primary is easy to track with the He II lines (right panel of Fig C.1), which are dominated by the hotter (A) component and to which the cooler (B) one contributes 10-15% of the total equivalent widths and is seen only as wings to the profile when the velocity separation is large. To see both components clearly we need to use one of the (weak) He I lines, where the much higher strength of the B component is sufficient to overcome the magnitude difference and have a similar (or even larger) equivalent width as the A component (left panel of Fig. C.1). However, there are few He I lines available and obtaining good-quality data is time expensive: each of the epochs in Fig. C.1 required 30-minute integrations on a 3.5 m telescope.

Fig. C.1. (left) Two LiLiMaRlin epochs of the He Iλ5875.6 and Na Iλλ5889.951,5895.924 region for the Bajamar star. A double Gaussian profile (green line) is fit to the He Iλ5875.6 absorption profile, with the blue vertical lines marking the A spectroscopic component and red vertical lines the B spectroscopic component. (right) The same two spectra in the region of the He IIλ8236.8 line. The blue and red vertical lines mark the expected positions as derived from the velocities of the He Iλ5875.6 fit. In both panels the spectral resolution has been degraded to R = 10 000 to better visualize the spectra. Telluric lines have been removed following Gardini et al. (2013).

We are currently obtaining GOSSS and LiLiMaRlin observations of the Bajamar star to derive its SB2 orbit but we still require further epochs before we can be certain about the result. Such an orbit will be included in one of the future installments of the MONOS project (Maíz Apellániz et al. 2019b; Trigueros Páez et al. 2021). Nevertheless, we can already present two preliminary results here: the period is around two weeks and the mass ratio is 2.2±0.5 (Fig. C.1).

Appendix D: TESS light-curve analysis

We have used TESS to analyze the photometric variability of the nine systems involved in ejection events from the Bermuda cluster. TESS light curves were extracted using the lightkurve (Cardoso et al. 2018) Python package (version 2.11), following the same procedure developed in Trigueros Páez et al. (2021). The data reveals only one eclipsing system, HD 200 776, plus five systems displaying pulsational or rotational variability (Bajamar, Toronto, Tyc 3179-00576-1, HD 195 965, and Tyc 3157-00918-1). The light curves for the other three systems, namely, Tyc 3575-01514-1 (observed in sectors 15 and 16), HDE 227 090 (sectors 14 and 15), and HD 201 795 (sector 15) display brightness stability at the σ = 0.5 mmag level for the first two stars and at the σ = 0.8 mmag for the last one.

HD 200 776 is an eclipsing binary (P  =  1.512 99 d) with an Algol-type light curve. Based on a spectroscopic and photometric analysis of the system, Dervişoǧlu et al. (2011) proposed that the system is a classical Algol system, where the secondary is filling the Roche lobe interacting with the primary component, that is, it is an evolved system. They derived absolute masses for the components of the system of M₁  =  6.1 M_⊙ and M₂  =  1.2 M_⊙ but they did so under the assumption of a distance of 501 ± 5 pc, which is excluded as a possibility by the Gaia EDR3 data. The source of the distance discrepancy is likely their assumption of a T_eff of 18 kK for the primary, which is too cool for a star of spectral type B1: the system must be more distant and massive. Their mass ratio and systemic radial velocity should be correct, however, as they are derived from the analysis of the radial velocity curve. The TESS light curve (Fig. D.1), observed in sectors 15 and 16, shows a similar morphology as the light curves presented by Dervişoǧlu et al. (2011), with a large difference in the depth between the primary and secondary eclipses and a continuous change in between them. We propose a different scenario to explain the shape of the light curve: a bloated secondary strongly irradiated by the hotter primary, a type of object that fits the measured mass ratio (Moe & Di Stefano 2015). This scenario makes the system compatible with an age of ∼2.0 Ma. In a forthcoming work, we will present a detailed modeling of the system.

Fig. D.1. Phased light curves for four of the objects analyzed with TESS. Symbols and colors are used to differentiate cycles. The number of epochs and cycles are given. The HJD of the zero point and the period can be read from the labels of the x axis.

The TESS observations of Tyc 3179-00756-1 correspond to sector 15. The light curve shows a small amplitude (∼1.5 mmag) regular signal with P  =  0.871 5 d (frequency, f = 1.147 42 d⁻¹). This type of coherent variability has been observed in main-sequence early-B type stars (Burssens et al. 2020) and it corresponds to slowly pulsating gravity modes with frequencies below 4 d⁻¹.

HD 195 965 was observed by TESS in sectors 15 and 16. The light curve (Fig. D.1) shows a strong and coherent periodic signal with P = 0.868 11 d (frequency, f = 1.151 93 d⁻¹), resembling a low-frequency pulsator (Burssens et al. 2020) but rotation modulation in a fast rotator cannot be discarded.

Tyc 3157-00918-1 was observed in sectors 14 and 15 and presents what is likely a rotational modulation with a period P = 6.928 d (Fig. D.1). Alternatively, if this object is a similar-mass early-A-type binary as suggested above, the value could correspond to one or one-half of the orbital period.

The Bajamar star was observed in sector 15. The TESS light curve shows stochastic low-frequency variability of hours scale superimposed on a small amplitude wave with P = 13.38 d. This value could correspond to the orbital period (see Appendix C) or it could be an artifact of the normalization procedure applied to the light curve.

TESS observed the Toronto star in sectors 15 and 16. The light curve (Fig. D.1) shows a very coherent signal with a period P = 4.682 d, with smaller amplitude stochastic variability. This photometric period is very different from the spectroscopic binary period of 48.5 d (Williams et al. 2001). The projected rotational velocity is 66 km/s (Simón-Díaz et al. 2017). Assuming a stellar radius of 7-8 R_⊙ we obtain a likely rotation period of about 4.6-5.3 d. Therefore, the photometric variability can be related to a rotation modulation of the O star.

All Tables

All Figures

Fig. 1. Top left: chart in equatorial coordinates of the three ejection events in this paper, color-coded in blue (Bajamar), red (Toronto), and dark green (HD 201 795). Colored solid lines show a representative trajectory for each system in the Sun’s LSR and colored short-dashed lines the equivalent after subtracting the motion of the Bermuda cluster. The black long-dashed line shows the Galactic equator and the square indicates the location of the bottom panels. The Cygnus asterism is also plotted for reference. Top right: DSS2 image of the same field as the top left panel. Trajectories are color-coded here in cyan (Bajamar), yellow (Toronto), and green (HD 201 795). Bottom left: equivalent chart of the square region in the top panels. Stars mark the location of each event in the Sun’s LSR (filled) and after subtracting the motion of the Bermuda cluster (non-filled). The ellipses mark the corresponding two-sigma uncertainty regions of each event using Gaussian weights. The black long-dashed line shows the Galactic equator and the black crosses mark the cluster members from Villafranca II. Bottom right: GALANTE (Maíz Apellániz et al. 2021a) three-color mosaic (red: F861M+F660N, green: F665N+F515N, blue: F450N+F420N) image of the same field as the bottom-left panel. The slight curvature of the trajectories is caused by the Galactic potential.

In the text

Fig. 2. Gaia EDR3 chart of the 144 members of the Bermuda cluster selected in Villafranca II. Symbol sizes encode G′ (corrected Gaia EDR3 magnitude), symbol colors encode $G_{\text{BP}}^{\prime }-G_{\text{RP}}^{\prime }$ (corrected Gaia EDR3 color), and arrows encode relative proper motions. Red arrows are used for the proper motions of the 17 brightest systems in G′. A cross indicates the used cluster center.

In the text

Fig. 3. Relative proper motions in right ascension (left) and declination (right) for the 144 members of the Bermuda cluster selected in Villafranca II as a function of the difference in each coordinate with respect to the cluster center. Red symbols are used for the 17 brightest systems in G′ and black symbols for the 127 fainter systems. The corresponding separation and velocity scales are given in the top and right axes, respectively. The correlation factors for the bright and faint systems are also given. The region in gray shows the possible linear fits to the faint systems in each case after rejecting outliers to leave 104 objects (left) and 107 (right).

In the text

Fig. A.1. Color-magnitude diagram with the extinction-corrected values (filled color circles) from Table A.2 and extinction-uncorrected ones (unfilled color circles). Some systems are labeled to allow for an easier correspondence with Table A.2. Different colors are used for the Bajamar (red), Toronto (blue), and HD 201 795 (dark green) events and for the eight brightest members of the Bermuda cluster (magenta). The Bajamar B and Toronto B components are not plotted but are discussed in the text. Unfilled black circles show the extinction-uncorrected measurements for the rest of the Bermuda cluster members selected in Villafranca II. The two black lines are the 1.78 Ma solar-metallicity Geneva isochrone with no rotation (a) without extinction and (b) with E(4405 − 5495) = 1.0 mag and R5495 = 3.0. Magnitudes are from Gaia EDR3.

In the text

	Fig. B.1. WISE W4+W3+W2 red-green-blue (RGB) mosaics for four walkaway and runaways systems. Each field is 20′×20′ (4.7 pc × 4.7 pc at a distance of 798 pc) and is oriented with north toward the top and east toward the left. In each mosaic the runaway candidate is at the center and the arrows show the absolute (green) and the relative-to-the cluster proper motions (yellow).
In the text

Fig. C.1. (left) Two LiLiMaRlin epochs of the He Iλ5875.6 and Na Iλλ5889.951,5895.924 region for the Bajamar star. A double Gaussian profile (green line) is fit to the He Iλ5875.6 absorption profile, with the blue vertical lines marking the A spectroscopic component and red vertical lines the B spectroscopic component. (right) The same two spectra in the region of the He IIλ8236.8 line. The blue and red vertical lines mark the expected positions as derived from the velocities of the He Iλ5875.6 fit. In both panels the spectral resolution has been degraded to R = 10 000 to better visualize the spectra. Telluric lines have been removed following Gardini et al. (2013).

In the text

	Fig. D.1. Phased light curves for four of the objects analyzed with TESS. Symbols and colors are used to differentiate cycles. The number of epochs and cycles are given. The HJD of the zero point and the period can be read from the labels of the x axis.
In the text