© ESO 2021

1. Introduction

Massive stars (M > 8 − 10 M_⊙) are hot and luminous and evolve rapidly (lifetimes up to a few tens of mega year). At the end of their lives, massive stars collapse producing a supernova and/or a gamma-ray burst. The end product is a compact object: a neutron star (NS) or a black hole (BH). During the supernova, various (heavy) chemical elements are produced; the outflow enriches the interstellar medium providing the building blocks for future generations of stars. For a review on the pre-supernova evolution of massive single and binary stars, see Langer (2012). For a review on the explosion mechanisms of core-collapse supernovae, we refer readers to Janka (2012).

The majority of massive stars are in binary (or multiple) systems (Sana et al. 2012). In a close binary, a phase of mass transfer before the supernova explosion of the primary results in an inversion of the mass ratio with the secondary becoming the most massive star in the system. This provides the condition that the binary system may remain bound after the supernova explosion; the latter also depends on the details of the natal kick (see, e.g., Kalogera et al. 1998). From that moment on, the binary consists of an OB-type star and a compact object that subsequently evolves into a high-mass X-ray binary (HMXB). X-rays are produced by the accretion of material from either the OB-star wind or through Roche-lobe overflow onto the compact object. HMXBs represent an important phase in the evolution of massive (close) binaries toward the formation of gravitational-wave sources: the merging event of a NS+NS, NS+BH or a BH+BH system (Belczynski et al. 2002). For a review on massive binary evolution, see Tauris & van den Heuvel (2006).

High-mass X-ray binaries are unique laboratories to test accretion physics and stellar evolution. HMXBs also provide the opportunity to accurately determine the physical properties of both the massive star and the compact object. If the compact object is an X-ray pulsar, its orbital parameters can be measured with high precision such that, for example, the masses of both companions are obtained given the radial-velocity amplitude of the OB star and the system inclination, such as in eclipsing systems. The X-ray eclipse duration provides a measure of the radius of the OB star. With the (Gaia) parallax, the OB-star spectrum and an estimate of the interstellar extinction, the effective temperature and luminosity of the OB star are obtained. As a consequence, the stellar parameters of the OB stars in HMXBs are amongst the most accurately known (cf. Kaper 2001).

Typically, due to their short lifetime, massive stars go into supernova close to the region where they were born. However, some massive stars move with high (supersonic) velocity through the interstellar medium, the so-called OB runaway stars (and walkaway stars) (Blaauw 1961; Renzo et al. 2019). They travel through interstellar space and when they explode as a supernova, they may have reached more remote gas-rich regions in our Galaxy where star formation can be triggered. Alternatively, they explode at a relatively high distance above the Galactic plane, such that the nuclearly enriched material can escape from the Galaxy.

In the vast majority of cases, the progenitor of the compact object was the originally most massive star in the binary (Pols 1994). This star depletes its hydrogen core faster than its companion, and expands into the supergiant phase first. If during this process the radius of the primary exceeds the Roche-lobe radius, it will transfer mass from its envelope onto the secondary. The secondary accretes hydrogen-rich material such that it becomes the most massive star in the system and the mixing of hydrogen into the core rejuvenates the star (van Bever & Vanbeveren 1998; Schneider et al. 2016; Hellings 1983). The primary becomes a helium star and produces a compact object after core collapse. If less than half of the total mass of the system is lost, the system can remain bound (Boersma 1961), depending on the natal kick, and becomes a HMXB when the accretion rate onto the compact object results in an observable X-ray flux.

Most HMXBs are runaway systems (van den Heuvel et al. 2000); when they move supersonically through the interstellar medium, a wind-bow shock can be produced (Kaper et al. 1997; Gvaramadze et al. 2011; Meyer et al. 2014; Prišegen 2019). Blaauw (1961) proposed that an OB runaway obtains its high space velocity as a consequence of a supernova explosion in a binary. If the system remains bound, the binary will obtain a recoil velocity comparable to the orbital velocity of the supernova progenitor. In this scenario, HMXBs are predicted to be runaway systems (unless no or little mass is lost from the system). However, the asymmetry of the supernova explosion likely introduces an additional velocity component (kick velocity) that may disrupt the binary, cf. Renzo et al. (2019). An alternative scenario to produce OB runaways is through the dynamical ejection from a stellar cluster (Poveda et al. 1967). This is more likely to occur in a young cluster when the stellar density is still high. Hoogerwerf et al. (2000) demonstrated that both mechanisms are at work; see also Jones et al. (2020) and Jilinski et al. (2010).

The time that passed since the OB-runaway system left its parent cluster is called the kinematical age of the system. If the kinematical age is known, it can be used to discriminate between the two scenarios. If an OB-runaway is produced by Blaauw’s scenario, the kinematical age and the (current) age of the parent cluster constrain the mass of the progenitor star that exploded during the supernova explosion (Ekström et al. 2012).

With a spectral type of O6.5 Iaf+, HD 153919 is the potentially most massive primary in the known sample of HMXBs in the Galaxy (cf. Falanga et al. 2015). The classification f+ refers to the presence of a number of emission lines in the spectrum (He II 4686, the N III complex at 4640 Å, and Si IV 4089,4116) produced by the dense stellar wind (Walborn 1990). The nature of the compact secondary has been subject of a long debate, e.g., Brown et al. (1996): it is either a massive neutron star (Clark et al. 2002) or a low-mass black hole, with a mass very similar to one of the merging compact objects in GW190412 (Abbott et al. 2020). No X-ray pulsations have been detected that would definitely identify 4U 1700-37 as a neutron star. An observational constraint on the mass ratio of the eclipsing system is obtained from the radial-velocity amplitude K and eccentricity of the orbit of HD 153919 (Hammerschlag-Hensberge et al. 2003).

Ankay et al. (2001) proposed that HD 153919/4U 1700-37 originates from the OB association Sco OB1 of which NGC 6231 is the suggested core (Fig. 1). Based on HIPPARCOS data they derive a kinematical age of 2 ± 0.5 million years. They take an upper limit on the age of Sco OB1 of 8 Myr and conclude that the progenitor of 4U 1700-37 went into supernova within 6 Myr, which puts a lower limit on the initial mass of the progenitor of 30 M_⊙.

Fig. 1. NGC 6231, the core ionizing cluster of the OB association Sco OB1, and the high-mass X-ray binary HD 153919/4U1700-37. The orange arrows show the distance the objects will travel in 0.5 Myr. The nearby open cluster Cl Trumpler 24, surrounded by the Prawn Nebula (IC 4621), and 4 bright stars in the Scorpius constellation are labeled. Source: Aladin.

NGC 6231 is a young open cluster that hosts many O and B stars and a large pre-main-sequence (PMS) population. It is considered to be the nucleus of the Sco OB1 association. The star formation process in this cluster has recently ended, and the molecular cloud in which the stars were formed has dispersed. This allows for high quality (optical) observations of the stars in the cluster and for a detailed study of the final outcome of the star formation process. The most recent age determinations of NGC 6231 are listed in Table 1.

NGC 6231 is the most probable birthplace of the HMXB 4U 1700-37. Feinstein et al. (2003) found evidence for a supernova explosion in NGC 6231 and showed that the observed initial mass function (IMF) allows for the initial presence of one or more very massive stars with masses up to about 70 M_⊙ (spectral type O3–O4). Sung et al. (2013) expect that NGC 6231 hosted 2-3 massive stars that were more massive than the currently most massive star in that cluster.

The astrometric and photometric data provided by the Gaia mission (Gaia Collaboration 2016, 2018, 2020) make it possible to quantitatively confirm the runaway scenario sketched by Ankay et al. (2001) based on HIPPARCOS data. The Gaia satellite executes an ambitious project of the European Space Agency (ESA) to record astrometric, photometric and spectroscopic data of more than one billion objects in the Galaxy and the Local Group up to G = 21 (Gaia Collaboration 2016; Riello et al. 2018). This is more than one percent of the stellar population in the Milky Way. With the accurate parallaxes and proper motions, we perform a membership analysis of NGC 6231. Subsequently, its current age is determined and the motion of HD 153919 is obtained relative to NGC 6231.

The scientific importance of confirming this scenario is that the progenitor mass of the compact object in 4U 1700-37 can be determined, the time of the supernova is set, the amount of material lost from the system is constrained, and the time it takes for the (rejuvenated) secondary, the current Of star HD153919, to evolve off the main sequence, can be measured (given the short duration of the HMXB phase). Thus, this provides a unique opportunity to test and constrain evolutionary models of massive binaries yielding binary neutron stars and/or black holes, the progenitors of the recently detected gravitational wave sources (e.g., Abbott et al. 2017).

In Sect. 2 we present the membership analysis of NGC 6231. We determine its age in Sect. 3. We provide an update of the parameters of the HMXB 4U1700-37 in Sect. 4. Section 5 presents the detailed reconstruction of the history of the system in space and time. In Sect. 6 we discuss different evolutionary scenarios that would fit these observations and address the physical nature of 4U1700-37.

2. Membership NGC 6231

The initial sample is manually selected from the Gaia DR2 database based on astrometric properties: coordinates, proper motion and parallax¹. We started with the candidate members of Sco OB1 listed in Ankay et al. (2001) to get a first impression of their astrometric properties. Subsequently, we queried the Gaia database at large in the thus defined parameter space. The resulting sample contains candidate members but also many (unrelated) stars in the field. After applying the astrometric corrections a membership analysis program is used to make a distinction between the field stars and the group members based on their proper motion and parallax Guo et al. (in prep.).

2.1. Manual selection

We start with a sample downloaded from the Gaia archive with an ADQL query (Appendix A). This is done by selecting a range of Galactic longitude, Galactic latitude and parallax based on our prior knowledge of the cluster. The dataset includes astrometry and photometry information of 127,349 sources. The next step is to further reduce the sample based on Galactic longitude, Galactic latitude, parallax and proper motion over-densities using visualization tools such as TOPCAT², until we obtain a sample that spans the cluster parameter space, but also includes enough field stars to make a distinction between cluster members and field stars. This is an important prerequisite for the membership analysis program. The remaining sample consists of 9101 sources. After astrometric corrections (Sect. 2.2.1 and Fig. 2), the remaining sample consists of 3661 sources. The positions and proper motions of the selected stars are shown in Figs. 6 and 7. A cut was made in distance such that only stars with a distance between 1300 pc and 1900 pc were included in the membership analysis.

Fig. 2. Flowchart showing the scripts that are used to get from the data to the results: yellow blocks refer to data; black blocks represent the applied manual selection criteria; green blocks indicate where astrometric corrections have been applied; red blocks represent the occurrence of photometric corrections; blue blocks signal the pipelines that are used to select the data and pink blocks represent the results.

2.2. Corrections to the Gaia DR2 data

In order to use the astrometric and photometric data provided by Gaia DR2 some corrections have to be made to the data obtained from the archive. Figure 2 shows all the necessary corrections that are needed to work with the data: astrometric corrections are displayed as green blocks while photometric corrections are shown as red blocks.

2.2.1. Astrometric corrections

After the release of Gaia DR2, errors on astrometric variables were found to be underestimated and some sources did not have reliable astrometric data at all. To correct for this all the errors on parallax and proper motion have to be recalculated and some sources have to be discarded from the sample.

External uncertainty correction. For the Tycho-Gaia Astrometric Solution (TGAS) of Gaia DR1 the published errors were obtained from a comparison with HIPPARCOS parallaxes. No such calibration has been done for Gaia DR2. This results in underestimated errors for all sources in Gaia DR2; faint sources (G > 13) being more affected than bright sources. Lindegren et al. (2018) tested the correctness of the errors provided by the Gaia catalog and came up with corrections on the errors of the parallaxes and the proper motions for bright sources (G < 13) and faint sources by comparing astrometric data from other catalogs, such as HIPPARCOS data and quasar data, to the astrometric data provided by Gaia DR2. The corrections that should be applied are shown in Eq. (1) in which σ_ext is the total error, σ_i is the error from the Gaia catalog and σ_s is a value computed by Lindegren et al. (2018) shown in Table 2. The value of k is 1.08 for all corrections.

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{ext}}=\sqrt{k^2{\sigma }_i^2+{\sigma }_s^2}.\end{array}$$(1)

The Renormalised Unit Weight Error (RUWE). This is a recommended goodness of fit indicator for Gaia DR2 astrometry that was not directly listed in the Gaia Archive³. It is a renormalization of the Unit Weight Error (UWE) by u₀(G, C): an empirical normalization factor dependent on photometric properties provided by ESA. If for a source RUWE > 1.4 holds, the goodness of fit is insufficient and it is advised to discard this source from the sample. The parameters needed to calculate UWE and tables with the correct normalization factor u₀(G, C) can be found on, respectively, the Gaia archive and the Gaia DR2 known issues web page⁴. The parameters that were used from the Gaia archive are the following:

• χ² = astrometric_chi2_al

• N = astrometric_n_good_obs_al

• G = phot_g_mean_mag

• C = bp_rp (if available)

UWE and RUWE are calculated using Eqs. (2) and (3):

$$\begin{array}{c}\hfill \mathrm{UWE}=\sqrt{{\chi }^2/(N-5)},\end{array}$$(2)

$$\begin{array}{c}\hfill \mathrm{RUWE}=\mathrm{UWE}/u_0(G,C).\end{array}$$(3)

No sources in our sample had to be discarded because they all have RUWE < 1.4.

Parallax and distance. The parallax is essential to determine the distance of the cluster as well as the runaway system, with which we are able to confirm whether the runaway system originates from the parent cluster in 3-dimensional space. However, despite Gaia’s unprecedented astrometric accuracy, the parallax measurements for these distant sources are relatively uncertain. Besides, the task of inferring distance from parallax with a properly defined uncertainty is not trivial and is dependent on the choice of a prior function on the distance distribution (Bailer-Jones et al. 2018). The naive method of simply inverting the parallax for distance is only acceptable for parallax_over_error > 5 (Bailer-Jones 2015). In this work, for simplicity, we use sources with parallax_over_error > 5 only, and invert them for the distance estimate (and subsequently updated to the EDR3 parallax).

2.2.2. Photometric corrections

The measured magnitudes in the GaiaG_BP and G_RP band may be contaminated such that these sources must be discarded from the sample. The G-band magnitude is corrected for all the remaining sources because of a systematic error in that band.

BP/RP excess filter. For the majority of the sources in the Gaia archive the flux is obtained in three bands: G, G_BP and G_RP. The flux in the G band is the most reliable because it is determined from profile-fitting to a narrow image. The G_BP and G_RP bands on the other hand give the total flux of the area of 7.35 arcsec² around the source which makes them more sensitive to pollution from the background or a nearby star. The G_BP and G_RP bands together cover about the same wavelength range and have similar transmission as the G band. The wavelength range that these passbands span and their transmission is shown in Fig. 3.

Fig. 3. The colored lines show the passbands for G (green), GBP (blue) and GRP (red), defining the Gaia DR2 photometric system. The thin, gray lines show the nominal, prelaunch passbands model published in Jordi et al. (2010), used for Gaia DR1. Source: ESA.

Evans et al. (2018) describe how to measure the pollution in the G_BP and G_RP bands and how to exclude bad measurements. First the assumption is made that if the G_BP and G_RP bands are not heavily polluted, the ratio of fluxes (I_GBP + I_GRP)/I_G, where I stands for flux, should be slightly bigger than one because the G_BP and G_RP bands have some overlap and a better transmission at some wavelengths, but together they cover almost the same area as the G band (Fig. 3). This ratio can be found in the Gaia archive as phot_bp_rp_excess_factor. If the phot_bp_rp_excess_factor significantly exceeds 1, then there must be pollution in the G_BP and/or G_RP band. A plot of (I_GBP + I_GRP)/I_G versus G_BP − G_RP clearly shows that some sources are polluted in the G_BP and G_RP band. Evans et al. (2018) computed a threshold that distinguishes the "well-behaved" sources from the polluted ones. If for a source: (I_GBP + I_GRP)/I_G < 1.3 + 0.06(G_BP−G_RP)² then it will be found under the red curve in Fig. 4 which means that the G_BP and G_RP band magnitude can be used. If it is situated above this curve then the G_BP and G_RP bands of this source are contaminated; this is the case for five sources in our sample.

Fig. 4. Flux excess versus color of sources for the members of NGC 6231 that were selected using our membership probability code. The red line corresponds to 1.3 + 0.06(GBP−GRP)2. The uncontaminated sources are below this curve (Evans et al. 2018). Five of the members of NGC 6231 are situated above the red curve; these sources were discarded from the sample for the age determination.

Correction of the G-band magnitude. The Gaia DR2 photometry in the G band is affected by systematic errors. It shows an approximately linear trend between G = 16 and G = 6 of 3.2 ± 0.3 mmag/mag (Maíz Apellániz & Weiler 2018). All the sources in the sample used for this study are in the range 6 < G < 16, and are corrected using:

$$\begin{array}{c}\hfill G^{\prime }=G-0.0032(G-6)\mathrm{for}6<G<16.\end{array}$$(4)

2.3. Membership probability

To distinguish cluster members from field stars a program is used to calculate the probability of a star being a member of a moving group based on a maximum-likelihood method which is a modified version of the method described by Lindegren et al. (2000). It makes use of a kinematic model.

2.3.1. The velocity model

First we assume that all stellar members of the cluster are moving in the same direction with the velocity vector:

$$\begin{array}{c}\hfill {\mathit{v}}_g={(v_{x,g},v_{y,g},v_{z,g})}^T.\end{array}$$(5)

We use ICRS coordinates in the calculation of the membership analysis. The next step is to use the normal triad $\widehat{p},\widehat{q},\widehat{r}$ to project the velocity vector v_g on the sky for a given sky coordinate (α, δ) and parallax ϖ to obtain a proper motion and radial velocity. The components of the normal triad are defined as:

$$\begin{array}{c}\hfill \widehat{p}=\left[\begin{array}{c}-sin\alpha \\ cos\alpha \\ 0\end{array}\right],\widehat{q}=\left[\begin{array}{c}-sin\delta cos\alpha \\ sin\delta sin\alpha \\ cos\delta \end{array}\right],\widehat{r}=\left[\begin{array}{c}cos\delta cos\alpha \\ cos\delta sin\alpha \\ sin\delta \end{array}\right].\end{array}$$(6)

Subsequently, we obtain the proper motion and radial velocity using the normal triad and the group motion:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mu }_{\alpha \ast ,g}\\ {\mu }_{\delta ,g}\\ v_{r,g}\\ \end{array}\right]=\left[\begin{array}{ccc}{p}_x\varpi /A& p_y\varpi /A& p_z\varpi /A\\ q_x\varpi /A& q_y\varpi /A& q_z\varpi /A\\ r_x& r_y& r_z\\ \end{array}\right]·\left[\begin{array}{c}{v}_{x,g}\\ v_{y,g}\\ v_{z,g}\end{array}\right],\end{array}$$(7)

where A  =  4.74047 yr km s⁻¹ is the astronomical unit in yr km s⁻¹, converting the angular unit of proper motion (mas yr⁻¹) to linear velocity (km s⁻¹). We assume a velocity dispersion model S that is the same across the whole group:

$$\begin{array}{c}\hfill S_g=\left[\begin{array}{ccc}{\sigma }_{x,g}^2& 0& 0\\ 0& {\sigma }_{y,g}^2& 0\\ 0& 0& {\sigma }_{z,g}^2\end{array}\right],\end{array}$$(8)

where the subscript g means group. Another assumption is that the dispersions in different directions are uncorrelated and thus we set all the terms except for the diagonal terms to zero. The field model is a copy of the group model. The only difference is the assumption that is made for the dispersion. We assume a small dispersion for the group model and a large dispersion for the field model.

2.3.2. The likelihood function

The next step is to apply the kinematic modeling method from Lindegren et al. (2000). The proper motion of each star i, μ_i = (μ_α * ,i, μ_δ, i)^T will be compared to both model predictions μ_i, g and μ_i, f: the group proper motion and the field proper motion on the position of the star, respectively, where we assume that the differences Δμ_i, k = μ_i − μ_k, in which k ∈ {g, f}, follows a Gaussian distribution. The likelihood function is:

$$\begin{array}{c}\hfill p_{i,k}({\mathit{\mu }}_k|{\mathit{\mu }}_i)=\frac{1}{2\pi |D_i|}exp(\mathrm{\Delta }{\mathit{\mu }}_{i,k}^T{D}_i^{-1}{\mathit{\mu }}_{i,k}/2),\end{array}$$(9)

in which D_i is the variance matrix:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_i=C_i& +\left[\begin{array}{cc}{\widehat{p}}_i^{T}S{\widehat{p}}_i& {\widehat{p}}_i^{T}S{\widehat{q}}_i\\ {\widehat{q}}_i^{T}S{\widehat{p}}_i& {\widehat{q}}_i^{T}S{\widehat{q}}_i\end{array}\right]{\left(\frac{{\varpi }_i}{\mathrm{A}}\right)}^2\hfill \\ \hfill & +\left[\begin{array}{cc}{({\widehat{p}}_i·{\mathit{v}}_{\mathit{k}})}^2& ({\widehat{p}}_i·{\mathit{v}}_{\mathit{k}})({\widehat{q}}_i·{\mathit{v}}_{\mathit{k}})\\ ({\widehat{p}}_i·{\mathit{v}}_{\mathit{k}})({\widehat{q}}_i·{\mathit{v}}_{\mathit{k}})& {({\widehat{q}}_i·{\mathit{v}}_{\mathit{k}})}^2\end{array}\right]{\left(\frac{{\sigma }_{\varpi ,i}}{\mathrm{A}}\right)}^2.\hfill \end{array}\end{array}$$(10)

It consists of three terms: the first term is the covariance matrix of the observed proper motion, the second term is the variance of proper motion caused by the velocity dispersion, and the third term is the variance caused by the uncertainty in parallax.

The total likelihood of the group and field model that is consistent with the proper motion of the ith star is:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_i={\lambda }_g·p_{i,g}({\mathit{\mu }}_g|{\mathit{\mu }}_i)+{\lambda }_f·p_{i,f}({\mathit{\mu }}_f|{\mathit{\mu }}_i).\end{array}$$(11)

The combination factors (λ_g, λ_f) satisfy λ_g + λ_f = 1. The program numerically maximizes the value of $\prod _{i=1}^N{\mathrm{\Phi }}_i$ for N sources by optimizing the parameters (v_g, S_g, v_f, S_f, λ_g, λ_f). When the optimization is finished, the probability that the ith star belongs to a group can be calculated by using the total likelihood as a normalization factor:

$$\begin{array}{c}\hfill p_{i,g}=\frac{{\lambda }_g·p_{i,g}({\mathit{\mu }}_g|{\mathit{\mu }}_i)}{{\mathrm{\Phi }}_i}.\end{array}$$(12)

Logically, the probability that the ith star belongs to the field is:

$$\begin{array}{c}\hfill p_{i,f}=\frac{{\lambda }_f·p_{i,f}({\mathit{\mu }}_f|{\mathit{\mu }}_i)}{{\mathrm{\Phi }}_i}=1-p_{i,g}.\end{array}$$(13)

Formally p_i, g > 0.5 is the lowest requirement for a candidate to be considered a member, however, one would use a threshold higher than 0.5 to achieve a clean sample in the color-magnitude diagram for isochrone fitting. In this work a source is considered a cluster member if p_i, g > 0.9.

2.4. The members of NGC 6231

273 stars in NGC 6231 are identified as members with a probability higher than 90%. The distribution of the membership probability of all candidates is shown in Fig. 5. Although a lower probability threshold would add more stars to the sample, it reduces the accuracy of the isochrone fitting. For more information on how the number of candidates changes with the probability threshold, see Table B.1 in Appendix B.

Fig. 5. Distribution of membership probability for the stars in NGC 6231 of all the candidates that were assessed by the membership analysis program. The first bin, containing over 2500 stars, is cut off at 100 stars. The bi-modal distribution shows that the program can distinguish members from field stars clearly, leaving only a small amount of ambiguous candidates in the middle range of probability. The black vertical line shows the threshold probability for membership.

Update to Gaia EDR3. In the light of Gaia’s Early Data Release 3 (Gaia EDR3, Gaia Collaboration 2020), which provides a significant improvement on the accuracy of parallax compared to Gaia DR2, we update the parallax and proper motions of the members in NGC 6231 and HD 153919/4U 1700-37 without rerunning the membership analysis based on Gaia DR2.

All the stars in the studied sample (magenta), and the selected members (blue) are shown in Fig. 6. The proper motions are displayed in Fig. 7. The clustering of the members in these figures clearly shows that the members are comoving in space.

Fig. 6. The studied sample projected on the sky (magenta points) and the sources that were selected as members (blue points).

Fig. 7. The studied sample in proper motion space (magenta points) and the sources that were selected as members (blue points). One can clearly see that the motion of the members is well confined. The proper motion of HD 153919 (green dot), a clear outlier from the cluster.

We estimate the distance by inverting the Gaia EDR3 parallax for the members in NGC 6231 as well as HD 153919, and plot the histogram in Fig. 8. The inverted mean parallax is at 1.63 kpc, in good correspondence with the distance of 1.64 and 1.59 kpc determined by Sana et al. (2006) and Sung et al. (2013), respectively. HD 153919/4U 1700-37 is also shown in Fig. 8 and appears to be at about the same distance as NGC 6231. In this work, we adopt 1.63 ± 0.15 kpc for the distance to NGC 6231.

Fig. 8. The distance of the members of NGC 6231 (blue histogram) and the most probable distance of HD 153919 (1.58 kpc, black line) with the corresponding distance interval (red and green lines). The average Gaia EDR3 distance to the cluster is 1.63 kpc, similar to the value determined by Sana et al. (2006) (1.64 kpc) and Sung et al. (2013) (1.59 kpc), both determined with the so-called spectroscopic parallax method (based on the spectral type of the star and the reddening).

3. The age of NGC 6231

With the cluster membership being established, we apply a photometric filter to the members. Five sources have too much pollution in the G_BP and G_RP bands and are thus discarded. Subsequently, the age determination is performed on 268 of the 273 members of NGC 6231. The determination of the cluster age was conducted by fitting isochrones to the members of the cluster in the (G − G_RP),G color-magnitude diagram. The method that was used is a simplified version of the isochrone fitting method by Jørgensen & Lindegren (2005). We use the code implemented by Guo et al. (in prep.).

3.1. Isochrone fitting

The age of NGC 6231 is determined with two different sets of isochrones (PARSEC and MIST), in order to check for consistency. The extinction A_V is treated as a free parameter. The Gaia catalog provides an estimate of A_G, but these values are not reliable and the conversion to the general extinction A_λ is not trivial (Jordi et al. 2010). The color excess E(B − V) varies with sight line, and thus the value for A_V. Sung et al. (2013) measured the total-to-selective extinction parameter R_V = 3.22 along this Galactic line of sight and E(B − V) varying between 0.45 and 0.60, with a slightly higher E(B − V) for some of the stars. The authors find an average E(B − V) of 0.47. This gives an A_V of 1.44 to 1.92 with an average of 1.50. Their field of view is somewhat different from ours. They observe a window of 0.67 deg by 0.67 deg while our field of view is a circle with a radius of 0.5 deg. Both sets of isochrones include 31 different values of A_V ranging from 0 to 4 mag in steps of 0.1. We adopt the solar metallicity: [Fe/H] = 0.

In this work we fit the isochrones to three sets of samples within the population of the cluster NGC 6231: the full sample, the pre-main-sequence (PMS) population (a subset of the full sample with (G − G_RP) > 0.65), and the main-sequence (MS) population (the complementary subset with (G − G_RP) < 0.65), so that the potential variation in age by analyzing different subpopulations is monitored. We do not take binarity into account.

3.1.1. PARSEC isochrones

For the PARSEC isochrones we apply a range from 0.2 Myr to 25 Myr with linear steps of 0.1 Myr. We used PARSEC release v1.2S + COLIBRI S_35. All the relevant isochrone parameters can be found in the Appendix in Sect. C.1 (Bressan et al. 2012; Chen et al. 2014, 2015; Marigo et al. 2017; Pastorelli et al. 2019; Tang et al. 2014).

3.1.2. MIST isochrones

For the MIST isochrones we use the range from 0.2 Myr to 25 Myr with steps of 0.1 Myr. We used MIST version 1.2. All the relevant isochrone parameters can be found in Appendix C.2 (Dotter 2016; Choi et al. 2016; Paxton et al. 2011, 2013, 2015).

3.2. Likelihood functions

On the color–absolute magnitude diagram (CMD) a 2-dimensional Gaussian function can be computed for every star with index i with color c_i, absolute magnitude m_i, and their errors σ_ci, σ_mi:

$$\begin{array}{c}\hfill g_i(c,m)=\frac{1}{2\pi {\sigma }_{c_i}{\sigma }_{m_i}}exp\{-\frac{1}{2}[{\left(\frac{c-c_i}{{\sigma }_{c_i}}\right)}^2+{\left(\frac{m-m_i}{{\sigma }_{m_i}}\right)}^2]\}.\end{array}$$(14)

Every point (c, m) on the CMD can be evaluated by g_i(c, m) and returns a probability of that point being the true color and magnitude of the star. Each isochrone of age a adjusted for extinction ϵ can be described as a curve I_a, ϵ(c, m). Now we can evaluate for every star i the probability for each isochrone I_a, ϵ(c, m) representing the age of the star

$$\begin{array}{c}\hfill G_i(a,ϵ)={\displaystyle \sum _j{g}_i(c_j,m_j),}\end{array}$$(15)

where j represents the index of the points on curve I_a, ϵ(c, m). Every star i now has a G-function G_i(a, ϵ) with a value for the likelihood for each isochrone I_a, ϵ(c, m). To calculate the age of a single stellar population, if we assume the formation process is coeval, we multiply the G-functions of all the stars to get one likelihood function with an evaluation of each isochrone used for the entire population:

$$\begin{array}{c}\hfill G(a,ϵ)={\displaystyle \prod _i{G}_i(a,ϵ).}\end{array}$$(16)

The maximum value of the likelihood function represents the best fit isochrone age and extinction for this cluster.

3.3. The age of NGC 6231

To determine the age of the stellar population of NGC 6231, the grid of PARSEC and MIST isochrones was used as described in Sects. 3.1 and 3.2.

PARSEC isochrone fitting. The results of the isochrone fitting with the grid of PARSEC isochrones are shown in Fig. 9. The plot shows a distribution based on the isochrone likelihood evaluation in the age-extinction space; the likelihood values are color coded (scale shown in the color bar). The location of the highest likelihood corresponds to the best fit isochrone of 4.7 Myr with an A_V of 1.7.

Fig. 9. G-function distribution in age and extinction based on the PARSEC isochrones for the whole sample. The likelihood is shown with a color scale indicated with the bar at the bottom of the figure. The best fit in this case is the isochrone of 4.7 Myr with an AV of 1.7 (see also Fig. 11).

MIST isochrone fitting. The results of the G-function of MIST isochrone fitting is shown in Fig. 10. The outcome is somewhat different from that obtained from the PARSEC isochrones: the best fit is the isochrone of 3.8 Myr with an A_V of 1.5, with a relatively large range in age.

Fig. 10. As Fig. 9: result obtained with the MIST isochrones for the whole sample. The best fit is the isochrone of 3.8 Myr with an AV of 1.5 (see also Fig. 11), with a relatively large range in age.

The best fitting isochrones based on both the PARSEC and MIST models are displayed in Fig. 11. Both models result in an acceptable fit to the data, however, with a slightly different outcome regarding the age. To investigate this further, we applied the same isochrone fitting procedure to subsamples of the population: (i) the PMS population; (ii) the MS population.

Fig. 11. The best fit MIST and PARSEC model isochrones to the full population, the main sequence only, and pre-main sequence only population, respectively. The absolute magnitude MG is based on the EDR3 parallax. It is clear that the isochrones that was fitted to the pre-main sequence only miss the fit due to the absence of the constraint from the MS stars, resulting in inconsistent extinction, regardless of models. Of the 4 remaining isochrones, the range in the age is between 3.8 Myr and 4.8 Myr and the AV that was determined is around 1.5–1.7 for all models, consistent with the literature values.

Age determination of the MS population. The best fit of the PARSEC isochrones to only the MS population of NGC 6231 results in an age of 4.8 Myr with an A_V of 1.7 mag. With the MIST isochrones we find an age of 4.4 Myr with an A_V of 1.5. Both results are similar to the fit on the whole population.

Age determination of the PMS population. Fitting the PARSEC and MIST isochrones to only the PMS population gives a similar age result of 4.5 and 4.8 Myr, respectively. However, both fits yield a higher extinction of 3.0 and 2.9 mag in A_V, double the literature value and our other results, with an A_V of 1.6. In this case, the best fitting isochrones miss the MS part of the population completely (Fig. 11) due to the lack of constraint from the MS sources in the CMD.

Estimating the uncertainties. We use a bootstrap method to estimate the uncertainty in the age. We repeat the age determination 200 times on randomly selected subsets of the sample. Both sets of isochrones are fit to a randomly selected sample that contains 90% stars of the original sample. For each sub-sample (All, PMS, MS), the bootstrap method returns the same result in the majority of the 200 trials. The estimated age and extinction of the whole and MS-only populations are consistent with only small standard deviation, while the PMS fittings show less stability with spurious results at very low age and extinction.

Table 3 and Fig. 12 show the results of the age and extinction determination with estimated uncertainty. In this work we adopt the results from the PARSEC isochrones fit to the whole sample, as the PARSEC model shows better consistency and stability in the results between the “all” and “MS” sample than the MIST model. In conclusion, our analysis results in an age and extinction of NGC 6231 of 4.7 ± 0.4 Myr and A_V = 1.7 ± 0.1 mag, respectively.

Fig. 12. Distribution of solutions based on the bootstrap method. The location of the circle represents the age and extinction of the best fit isochrone, and the size of the cirle is proportional to the occurrence rate based on the 200 trials.

4. Updating the system parameters of HD 153919/4U 1700-37

The new distance determination of HD 153919, the optical counterpart of the high-mass X-ray binary 4U 1700-37, provides the opportunity to reassess the system parameters. Based on the Gaia EDR3 parallax, the distance to the system becomes 1.58 ± 0.07 kpc (we note that the DR2 parallax resulted in 1.75 ± 0.24 kpc). The effective temperature corresponding to the O6.5 Iaf+ spectral type is 35, 000 ± 1000 K (Clark et al. 2002). With the X-ray eclipse duration the radius of the O supergiant can be determined: R_p = 25.1 ± 4.0 R_⊙, using an inclination i = 62 deg ± 1 deg and an eclipse semiangle θ_E = 32 deg ± 1 deg (Falanga et al. 2015). The (black-body) luminosity of HD 153919 becomes log(L/L_⊙) = 5.93 ± 0.05. This is consistent with the (astrometric) luminosity of HD153919, taking into account the photometry, bolometric correction, and extinction.

Clark et al. (2002) used Monte Carlo simulations to determine the mass of HD 153919 and 4U 1700-37: M₁ = 58 ± 11.5 M_⊙ and M_X = 2.44 ± 0.27 M_⊙, respectively. The latter mass is relatively high for a neutron star, but low for a stellar-mass black hole. Clark et al. (2002) report that the Monte Carlo simulations fit a zero eccentricity solution and a solution with an eccentricity of 0.22 equally well. The zero eccentricity solution is the conservative approach, because the masses of the O supergiant and compact object would be significantly larger with an eccentric orbit. They state that masses that high do not match the spectral type and log g of the O supergiant.

Hammerschlag-Hensberge et al. (2003) favor an eccentric orbit of e = 0.22. They state that Clark et al. (2002) miscalculated the masses of the binary members because the value for K (20.6 km s⁻¹) that was used by them was determined using an eccentric orbit; for a circular orbit they should have used K = 18.7 km s⁻¹. If they would have used the correct value of K for the circular orbit, then the masses would have been even higher than those proposed by Clark et al. (2002) for the eccentric solution. Hammerschlag-Hensberge et al. (2003) also found a trend in the residuals of K while trying to fit a circular orbit that could be solved when eccentricity was included. Hammerschlag-Hensberge et al. (2003) determine that the masses of the members of the binary are about 4% higher than determined by Clark et al. (2002) with K = 20.6 and e = 0.22. They calculated for HD 153919/4U 1700-37 M_P = 60 ± 11 and M_X = 2.54 ± 0.27. With the masses and mass function obtained by Hammerschlag-Hensberge et al. (2003) we find sin(i) = 0.88 (±0.14), i = 62. This corresponds well with the value that was determined by Falanga et al. (2015) (i = 62 deg ± 1 deg). The latter authors arrive at somewhat lower masses, but we adopt the masses of Hammerschlag-Hensberge et al. (2003).

A list of the updated parameters is given in Table 4. These parameters are internally consistent and will be used during the remainder of the paper.

5. The kinematical age of the runaway system

The aim of this section is to confirm that Sco OB1, and more specifically NGC 6231, is the parent OB association of the HMXB HD 153919/4U 1700-37 by reconstructing the path of the system based on Gaia DR2 data. In this way, the kinematical age of the system, and thus the time of supernova, can be determined. All members of NGC 6231 as well as HD 153919 have well measured proper motions. The radial velocities of the stars in NGC 6231 are not available in Gaia DR2, additionally, Sartoretti et al. (2018) state that the radial velocities of stars in DR2 with T_eff > 7000 K are unreliable. Therefore, radial velocities from earlier studies are used in the reconstruction of the path. For simplicity, the Galactic potential has not been taken into account when calculating the path.

5.1. Reconstruction of the path

To trace back the paths of the individual stars the coordinates, parallaxes, and proper motions were converted into a Cartesian coordinate system. The transformation of the coordinates is done with the following equations:

$$\begin{array}{c}\hfill x_0=dcosbcosl\end{array}$$(17)

$$\begin{array}{c}\hfill y_0=dcosbsinl\end{array}$$(18)

$$\begin{array}{c}\hfill z_0=dsinb,\end{array}$$(19)

where d is the distance, l is the Galactic longitude and b is the Galactic latitude. Section 2.3.1 describes how to transform vectors from a spherical coordinate system to a Cartesian coordinate system using the normal triad $(\widehat{p},\widehat{q},\widehat{r})$, see Eqs. (6) and (7). To transform vectors from a spherical coordinate system to a Cartesian coordinate system we have to multiply the proper motions and radial velocity with the inverse of the normal triad matrix ${(\widehat{p},\widehat{q},\widehat{r})}^{-1}$. The inverse of matrices representing a spatial rotation operation is the same as the transposed matrix ${(\widehat{p},\widehat{q},\widehat{r})}^T$. The resulting operation is the following:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]=\left[\begin{array}{ccc}{p}_{x}dA& q_{x}dA& r_x\\ p_{y}dA& q_{y}dA& r_y\\ p_{z}dA& q_{z}dA& r_z\\ \end{array}\right]·\left[\begin{array}{c}{\mu }_{\alpha \ast ,g}\\ {\mu }_{\delta ,g}\\ v_{r,g}\\ \end{array}\right],\end{array}$$(20)

with d the distance and A = 4.74047 yr km s⁻¹ to convert from km s⁻¹ to AU year⁻¹. After transforming the coordinates and velocities from spherical coordinates to Cartesian coordinates it is straightforward to trace the positions of the sources (x, y, z) back in time t.

$$\begin{array}{c}\hfill x=x_0+v_{x}t\end{array}$$(21)

$$\begin{array}{c}\hfill y=y_0+v_{y}t\end{array}$$(22)

$$\begin{array}{c}\hfill z=z_0+v_{z}t\end{array}$$(23)

where (x₀, y₀, z₀) is the position of the sources today, and the time t < 0 as we move the positions back in time. After obtaining the traced back Cartesian coordinates of the sources one can easily transform them back into the Galactic coordinate system:

$$\begin{array}{c}\hfill l={tan}^{-1}(y/x)\end{array}$$(24)

$$\begin{array}{c}\hfill b={tan}^{-1}(z/\sqrt{x^2+y^2})\end{array}$$(25)

$$\begin{array}{c}\hfill d=\sqrt{x^2+y^2+z^2}\end{array}$$(26)

5.2. Radial velocities

5.2.1. The radial velocity of NGC 6231

We obtained the radial velocity of seven OB stars in NGC 6231 from the literature (see Table 5) and used the mean value to represent the radial velocity of the cluster. All these OB stars have a binary companion except for HD 152314; HD 152270 is a Wolf-Rayet + O star binary. The listed radial velocities of the binaries are based on multiple measurements at different orbital phases of the binaries and the orbital motion is removed. The mean radial velocity of the cluster is −33.5 ± 4.0 km s⁻¹, in agreement with the measurement of −27.3 km s⁻¹ by Kharchenko et al. (2005).

5.2.2. Radial velocity of HD 153919

The determination of the radial velocity of HD153919 is complicated as this Of supergiant has a dense, supersonically expanding stellar wind and is in a binary system. Hutchings (1974) measured the velocity of multiple spectral lines and found it to be −60 km s⁻¹. We adopt the value of −64.5 ± 1.5 km s⁻¹ from Gies (1987), see also Kaper et al. (1994). The radial-velocity amplitude due to the binary motion is 20.6 ± 1.0 km s⁻¹ (Hammerschlag-Hensberge et al. 2003).

5.2.3. The space velocity of HD 153919

The relative velocity of HD 153919 with respect to the parent cluster NGC 6231 thus becomes v_sys = 63 ± 5 km s⁻¹. In the supernova scenario the runaway velocity provides a constraint on the amount of mass lost from the system during the supernova.

5.3. The kinematical age of 4U 1700-37

The path of HD 153919/4U 1700-37 relative to NGC 6231 has been reconstructed to determine the location of the closest encounter and the time that has passed since then: the kinematical age. The coordinates, proper motions, distances, the average radial velocity of 7 likely members of NGC 6231, and the radial velocity of HD 153919/4U 1700-37 (Table 5) were used to project the path on different sky planes: in Galactic coordinates l and b (Fig. 13) and in Cartesian coordinates x, y, and z (Fig. 14).

Fig. 13. Location of HD 153919/4U 1700-37 and the individual members of NGC 6231 now (blue) and 2.2 Myr ago (green). The orange arrows indicate the movement of HD 153919 and NGC 6231 in 0.5 Myr, respectively.

Fig. 14. The position of HD 153919 and that of the individual members of NGC 6231 now (blue) and 2.2 Myr ago (green) in Galactic Cartesian coordinates. The uncertainty in the distance of HD153919 is represented by the magenta lines. The orange arrows show the movement of HD 153919 and NGC 6231 in 0.5 Myr.

It becomes clear from Fig. 13 that, projected on the sky, HD 153919/4U 1700-37 was positioned in the outskirts of NGC 6231 2.2 ± 0.1 Myr ago. Figure 14 shows that, in the radial direction, the uncertainty in the distance of HD153919 (magenta error bar) is large. This is also reflected by the elongated shape of the cluster in the radial direction.

6. Discussion

From the previous sections we conclude that 2.2 Myr ago the HMXB HD 153919/4U 1700-37 escaped with a space velocity of 63 km s⁻¹ from the outskirts of NGC 6231, the 4.7 Myr young cluster in the center of the OB association Sco OB1. In this section we discuss the implications of the reported astrometric solution, propose scenarios to explain the evolutionary history of the system, and consider the physical nature of the compact object.

6.1. The astrometric solution

Gaia determines accurate parallaxes and proper motions, but the errors on the astrometric parameters increase with distance, and the distance to HD 153919/4U 1700-37 and the cluster NGC 6231 is relatively large (1.6 kpc). Figure 14 shows that many members of NGC 6231 appear to be up to 0.3 kpc away from the center of NGC 6231 in the radial direction. The elongated shape that the cluster seems to have is a consequence of the uncertainty in the distance determination for the individual members. If we assume the cluster to be spherical, we arrive at a physical diameter of about 20 pc.

In Fig. 13 the cluster seems to be more compact 2.2 Myr ago than today. To investigate whether this is the result of a physical rather than an apparent expansion, we calculated the average angular distance of the members to the center of the cluster in degrees during the past 5 Myr in steps of 0.1 Myr and corrected for the projection effect by the varying distance (the cluster is approaching us). The minimum average angular distance of the cluster members to the center is reached 1.8 Myr ago (0.19 deg), but the difference to the current average distance (0.21 deg) is small. Therefore, we cannot conclude that NGC 6231 has significantly expanded since its formation about 4.7 Myr ago.

One may naively expect that HD 153919/4U 1700-37 would originate from the center of NGC 6231 since the progenitors of 4U 1700-37 and HD 153919 were likely both (very) massive stars (cf. Sect. 6.2). The reconstruction of its path indicates that the binary was ejected from the outskirts of the cluster, which would not necessarily be the region with the highest stellar density for which the dynamical ejection scenario is the most efficient (e.g., Fujii & Portegies Zwart 2011; Oh & Kroupa 2016).

We conclude that HD 153919/4U 1700-37 is a runaway system that can be traced back to NGC 6231, its parent cluster. The current membership and age of NGC 6231, still including a number of massive stars, allows some stars to have been massive enough to have already ended their evolution (Sung et al. 2013; Feinstein et al. 2003). There is evidence for the occurrence of a supernova explosion in NGC 6231 (Feinstein et al. 2003), but it is unclear whether this event has a relation with the formation of 4U 1700-37.

6.2. Possible progenitor systems

We can use the parameters of HD 153919/4U 1700-37 to constrain some of the characteristics of the progenitor system. Nelemans et al. (1999) derived a relation between the mass lost during a supernova explosion (ΔM) and the space velocity v_sys (assuming only a Blaauw kick):

$$\begin{array}{c}\hfill \left(\frac{\mathrm{\Delta }M}{M_{\odot }}\right)=\left(\frac{v_{\mathrm{sys}}}{213{\mathrm{km}\mathrm{s}}^{-1}}\right){\left(\frac{M}{M_{\odot }}\right)}^{-1}{\left(\frac{P_{\mathrm{cir}}}{\mathrm{day}}\right)}^{1/3}{\left(\frac{M+m}{M_{\odot }}\right)}^{5/3},\end{array}$$(27)

where M is the present mass of HD 153919: 60 M_⊙, m is the mass of the compact object in 4U 1700-37: 2.54 M_⊙, and P_cir is the orbital period after recircularization of the orbit. We neglected the eccentricity to be able to use this equation with P_cir of 3.41 days. With v_sys of 63 km s⁻¹, this yields a ΔM of 7.4 M_⊙. This means that the core of the progenitor of 4U1700-37 ($M_1^{\text{pre-CC}}$) was about 7.5 M_⊙ after the RLOF phase and before the supernova. Such a core would be too small for a primary star with M > 60 M_⊙ (e.g., Brott et al. 2011), suggesting that the companion star HD153919 has accreted mass before the end of the main sequence of the progenitor of the compact object (e.g., Belczynski & Taam 2008), that is case A mass transfer (Kippenhahn & Weigert 1967).

Renzo et al. (2019) have described an analytical method to probe the change in orbital separation using the masses of the primary and secondary before and after Roche-lobe overflow. This was originally designed to treat case B RLOF, i.e. mass transfer after a He core is well established in the donor. While this might not be the case for 4U 1700-37, we proceed to demonstrate that assuming this scenario leads to inconsistent results.

We start with the assumption that the primary loses its complete envelope after mass transfer:

$$\begin{array}{c}\hfill M_1^{\mathrm{pre}-\mathrm{CC}}\simeq M_1^{\mathrm{ZAMS}}-{\mu }_{\mathrm{env}}{M}_1^{\mathrm{ZAMS}},\end{array}$$(28)

where $M_1^{\text{ZAMS}}$ is the initial mass of the primary and μ_env is the fraction of the primary that is the envelope. We rewrite to get:

$$\begin{array}{c}\hfill {\mu }_{\mathrm{env}}\simeq \frac{M_1^{\mathrm{ZAMS}}-M_1^{\mathrm{pre}-\mathrm{CC}}}{M_1^{\mathrm{ZAMS}}}.\end{array}$$(29)

Then for the secondary we have:

$$\begin{array}{c}\hfill M_2^{\mathrm{pre}-\mathrm{CC}}\simeq M_2^{\mathrm{ZAMS}}+{\beta }_{\mathrm{RLOF}}{\mu }_{\mathrm{env}}{M}_1^{\mathrm{ZAMS}},\end{array}$$(30)

where β_RLOF is the fraction of the envelope of the primary that is accreted by the secondary. We can rewrite this to:

$$\begin{array}{c}\hfill {\beta }_{\mathrm{RLOF}}\simeq \frac{M_2^{\mathrm{pre}-\mathrm{CC}}-M_2^{\mathrm{ZAMS}}}{{\mu }_{\mathrm{env}}{M}_1^{\mathrm{ZAMS}}}.\end{array}$$(31)

Since we know that $M_1^{\text{pre-CC}}$ ≃ 10 M_⊙ and $M_2^{\text{pre-CC}}$ ≃ 60 M_⊙, filling in $M_1^{\text{ZAMS}}$ and $M_2^{\text{ZAMS}}$ gives us values of μ_env and β_RLOF. If we find values that are physically comfortable we can proceed to calculate the orbital separation of the progenitor system. Assuming a symmetrical supernova that does not change the orbital separation we assume a_{pre − CC} to be the current orbital separation of 38 R_⊙. For the orbital separation of the progenitor system (a_ZAMS) we can write:

$$\begin{array}{c}\hfill a_{\mathrm{ZAMS}}=a_{\mathrm{pre}-\mathrm{CC}}(\frac{M_1^{\mathrm{pre}-\mathrm{CC}}{M}_2^{\mathrm{pre}-\mathrm{CC}}}{M_1^{\mathrm{ZAMS}}{M}_2^{\mathrm{ZAMS}}}{)}^2(\frac{M_1^{\mathrm{ZAMS}}+M_2^{\mathrm{ZAMS}}}{M_1^{\mathrm{pre}-\mathrm{CC}}+M_2^{\mathrm{pre}-\mathrm{CC}}}{)}^{2{\gamma }_{\mathrm{RLOF}}+1},\end{array}$$(32)

where γ_RLOF is the specific angular momentum of the matter that leaves the system. For simplicity, we assume γ_RLOF = 1, i.e. mass not accreted during the RLOF takes away the specific orbital angular momentum of the binary (Dominik et al. 2012). In reality, this parameter could also vary during the mass transfer phase.

To find out whether we can produce a system like HD 153919/4U 1700-37, we tested values of $M_1^{\text{ZAMS}}$ ranging from 30 to 80 M_⊙ and $M_2^{\text{ZAMS}}$ ranging from 10 to 50 M_⊙, with a γ_RLOF of 1, 2 and 3. A γ_RLOF of 2 or 3 seems unlikely, since this would mean that the system suffers extreme angular momentum losses. The results of these tests are displayed in Table 6. Only solutions with a high γ_RLOF and a low β_RLOF, so when a lot of angular momentum leaves the system due to nonconservative mass transfer, seem to give a a_ZAMS that could be large enough to avoid merging of the system. For comparison: HD 153919 with a mass of 60 M_⊙ has a radius of 22 R_⊙. This result is at least in part due to our assumption of taking the present-day separation as pre-CC separation and assume that the explosion does not change it.

6.2.1. Binary evolution channels

With the orbital parameters of HD 153919/4U 1700-37, the current mass of HD 153919 and 4U 1700-37, the space velocity of HD 153919/4U1700-37 relative to NGC 6231 and the determined lower limit of the mass of the progenitor of 4U 1700-37, population synthesis can be used to constrain the parameters of the progenitor system. The analytical approach in Sect. 6.2 showed that constraining the parameters of the progenitor system is not an easy task, as also pointed out before by Ankay et al. (2001).

If we assume conservative mass transfer, the effect of the mass transfer on the orbital separation can be approximated by:

$$\begin{array}{c}\hfill \frac{\dot{a}}{a}=2(\frac{M_d}{M_a}-1)\frac{\dot{M_d}}{M_d}.\end{array}$$(33)

This tells us that because Ṁ_d < 0 the orbit shrinks (ȧ < 0) as long as M_d > M_a but when the mass ratio flips and M_d < M_a, the orbit will expand (ȧ > 0) a lot more than it shrank before the mass ratio flipped.

Case A mass transfer. If we assume case A mass transfer, which overall leads to more conservative mass transfer because the mass transfer rate is relatively low and proceeds on a nuclear timescale, the orbit widens a lot as described above and the binary ends with a large orbital separation which does not fit our scenario well (assuming a_{pre − CC} = a_today).

Case B mass transfer. If we assume case B mass transfer, which happens when the donor has a well-established He core and fills its Roche lobe because it expands significantly, the mass transfer rate is way higher than with case A mass transfer. This generally means that mass transfer is less conservative. This leads to a more compact binary, but now a different problem emerges. A supergiant, by definition, depleted the hydrogen in its core and thus has a heavy helium core after the mass transfer phase, which would probably produce a black hole after the core-collapse that is way more massive than what we have observe for 4U 1700-37 (2.54 M_⊙).

A possible scenario. With some fine tuning one can create systems similar to what we observe. The most probable scenario includes case A mass transfer. An asymmetric supernova explosion with a very specific kick direction and amplitude can solve the complication of the widening of the binary (Wongwathanarat et al. 2013; Janka 2013, 2017). If the kick is directed outward and kicks the compact object toward its companion, then the orbital separation can shrink enough to reproduce the observed system.

6.2.2. Asymmetric supernova kicks

The discovery of radio pulsars with an extremely high space velocity exceeding in some cases 1000 km s⁻¹ (e.g Cordes et al. 1993) is best explained with an asymmetric supernova kick. The mechanism behind supernova kicks is not fully understood. If a supernova fails, most or all matter will fall back onto the core of the star and a black hole forms. It is thought that this happens if a very massive single star undergoes a supernova explosion and the kick received here is thought to be non to little since no matter is ejected and no extra velocity is needed to conserve momentum. Black-hole systems like Cyg X-1 (Mirabel & Rodrigues 2003) have low space velocities, which can be understood in terms of a failed supernova leading to fall back and the formation of a black hole, and no mass loss from the system. However, the low-mass black-hole system 4U1957+11 may have a high space velocity (Maccarone et al. 2020).

In order for a supernova explosion not to fail, it is thought that asymmetry plays a major role. The stellar plasma does not propagate outwards as fast as with a symmetric outflow because it is not accelerated only in the radial direction. This causes the density of the material to be high enough to absorb a big part of the energy carried away by the neutrino flow and accelerate outwards.

Janka (2013) explains the supernova kick with a so-called tug-boat mechanism in which the most massive part of the supernova ejecta moves out slower than the rest of the ejecta because of inertia, and pulls the neutron star gravitationally giving it a velocity in that direction. The faster outflow, that causes most of the explosive nucleosynthesis, would then be situated at the other side of the core. Whether (low mass) BHs can receive kicks remains an open question (Repetto et al. 2012; Mandel 2016; Janka 2017; Atri et al. 2019).

6.3. The nature of 4U 1700-37

It is not clear whether 4U 1700-37 is a neutron star or a black hole. It emits in X-rays but does not show significant X-ray pulsations identifying it as a neutron star. White et al. (1983) states that it is possible that the neutron star undergoes spherical accretion instead of accretion along the magnetic field lines. This could be related to the accretion mechanism. Wind accretion, rather than Roche-lobe overflow, may not lead to the formation of an accretion disk and thus give rise to spherical accretion such that X-ray pulsations are not observed. Alternatively, the spin and magnetic axis of the neutron star could be aligned.

The mass of 4U 1700-37 that is estimated by Hammerschlag-Hensberge et al. (2003) of 2.54 ± 0.27 M_⊙ is at the high end of the mass range observed for neutron stars (2 M_⊙), but lower than than that of the lowest mass black holes (4–5 M_⊙) observed in an X-ray binary system (e.g., Orosz et al. 2003). The X-ray spectrum of 4U 1700-37 shows characteristics of a neutron star (Seifina et al. 2016; White et al. 1983). Interestingly, LIGO/Virgo have recently announced the discovery of GW190412 (Abbott et al. 2020) which involved a compact object of similar mass and unknown nature. While the metallicity of 4U1700-37 is probably different from the metallicity of the progenitor of GW190412, we suggest that 4U1700-37 is a possible analog of its progenitor. Thus studying this system might shed light on the evolution of GW progenitors.

Perhaps the strongest argument in favor of a neutron star is that the observed runaway velocity indicates that a lot of material is lost from the system during the supernova explosion. However, it is unclear whether a low mass BH might receive a similar kick amplitude at birth.

The high space velocity of the system and the dense stellar wind are expected to result in the formation of a wind bow shock such as observed in Vela X-1 (Kaper et al. 1997) and 4U1907+09 Gvaramadze et al. (2011). Such a wind bow shock is not observed for 4U1700-37. However, only a fraction of the OB runaway stars produces a wind bow shock (Huthoff & Kaper 2002). This may be either due to the low density or high temperature of the ambient medium, or both.

Acknowledgments

We acknowledge an anonymous referee who provided constructive remarks that helped to improve the quality of the paper. We want to thank Jari van Opijnen, Mark Snelders, Rob Farmer and Selma de Mink for interesting discussions. This research has been supported by NOVA. 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.

References

Appendix A: ADQL selection candidate members NGC 6231

SELECT designation, source_id,
ra, ra_error, dec, dec_error, parallax,
parallax_error, parallax_over_error,
pmra, pmra_error, pmdec, pmdec_error,
ra_dec_corr, ra_parallax_corr, ra_pmra_corr,
ra_pmdec_corr, dec_parallax_corr, dec_pmra_corr,
sdec_pmdec_corr, parallax_pmra_corr,
parallax_pmdec_corr, pmra_pmdec_corr,
r_est, r_lo, r_hi,
astrometric_primary_flag, duplicated_source,
radial_velocity, radial_velocity_error,
phot_g_mean_flux, phot_g_mean_flux_error,
phot_g_mean_mag,
phot_bp_mean_flux, phot_bp_mean_flux_error,
phot_bp_mean_mag,
phot_rp_mean_flux, phot_rp_mean_flux_error,
phot_rp_mean_mag,
phot_bp_rp_excess_factor, bp_rp, bp_g, g_rp,
phot_variable_flag,
l, b, teff_val,
teff_percentile_lower, teff_percentile_upper,
a_g_val, a_g_percentile_lower,
a_g_percentile_upper,
e_bp_min_rp_val,
e_bp_min_rp_percentile_lower,
e_bp_min_rp_percentile_upper,
radius_val, radius_percentile_lower,
radius_percentile_upper,
lum_val, lum_percentile_lower,
lum_percentile_upper,
astrometric_chi2_al, astrometric_n_good_obs_al

FROM external.gaiadr2_geometric_distance
JOIN gaiadr2.gaia_source USING (source_id)
WHERE l >= 342.5 AND l <= 344
AND b >= 0 AND b <= 2 AND
parallax >= 0.4 AND parallax <= 1

Appendix B: Membership probability thresholds for NGC 6231

Appendix C: Isochrone parameters

C.1. PARSEC isochrone parameters

Thermal pulse cycles included

On RGB, assumed Reimers mass loss with efficiency eta = 0.2

Photometric system: Gaia’s DR2 G, G_BP and G_RP (Vegamags, Gaia passbands from Maiz-Apellaniz and Weiler 2018)

Using YBC version of bolometric corrections as in Chen et al. (2019)

O-rich circumstellar dust ignored

C-rich circumstellar dust ignored

IMF: chabrier_lognormal_salpeter

Kind of output: isochrone tables

C.2. Mist isochrone parameters

MIST version number = 1.2

MESA revision number = 7503

photometric system = UBV(RI)c, 2MASS, Kepler, Hipparcos, Gaia (Vega)

Yinit = 0.2703 Zinit = 1.42000E-02 [Fe/H] = 0 [a/Fe] = 0 v/vcrit = 0

All Tables

All Figures

	Fig. 1. NGC 6231, the core ionizing cluster of the OB association Sco OB1, and the high-mass X-ray binary HD 153919/4U1700-37. The orange arrows show the distance the objects will travel in 0.5 Myr. The nearby open cluster Cl Trumpler 24, surrounded by the Prawn Nebula (IC 4621), and 4 bright stars in the Scorpius constellation are labeled. Source: Aladin.
In the text

Fig. 2. Flowchart showing the scripts that are used to get from the data to the results: yellow blocks refer to data; black blocks represent the applied manual selection criteria; green blocks indicate where astrometric corrections have been applied; red blocks represent the occurrence of photometric corrections; blue blocks signal the pipelines that are used to select the data and pink blocks represent the results.

In the text

	Fig. 3. The colored lines show the passbands for G (green), GBP (blue) and GRP (red), defining the Gaia DR2 photometric system. The thin, gray lines show the nominal, prelaunch passbands model published in Jordi et al. (2010), used for Gaia DR1. Source: ESA.
In the text

	Fig. 4. Flux excess versus color of sources for the members of NGC 6231 that were selected using our membership probability code. The red line corresponds to 1.3 + 0.06(GBP−GRP)2. The uncontaminated sources are below this curve (Evans et al. 2018). Five of the members of NGC 6231 are situated above the red curve; these sources were discarded from the sample for the age determination.
In the text

Fig. 5. Distribution of membership probability for the stars in NGC 6231 of all the candidates that were assessed by the membership analysis program. The first bin, containing over 2500 stars, is cut off at 100 stars. The bi-modal distribution shows that the program can distinguish members from field stars clearly, leaving only a small amount of ambiguous candidates in the middle range of probability. The black vertical line shows the threshold probability for membership.

In the text

	Fig. 6. The studied sample projected on the sky (magenta points) and the sources that were selected as members (blue points).
In the text

	Fig. 7. The studied sample in proper motion space (magenta points) and the sources that were selected as members (blue points). One can clearly see that the motion of the members is well confined. The proper motion of HD 153919 (green dot), a clear outlier from the cluster.
In the text

Fig. 8. The distance of the members of NGC 6231 (blue histogram) and the most probable distance of HD 153919 (1.58 kpc, black line) with the corresponding distance interval (red and green lines). The average Gaia EDR3 distance to the cluster is 1.63 kpc, similar to the value determined by Sana et al. (2006) (1.64 kpc) and Sung et al. (2013) (1.59 kpc), both determined with the so-called spectroscopic parallax method (based on the spectral type of the star and the reddening).

In the text

	Fig. 9. G-function distribution in age and extinction based on the PARSEC isochrones for the whole sample. The likelihood is shown with a color scale indicated with the bar at the bottom of the figure. The best fit in this case is the isochrone of 4.7 Myr with an AV of 1.7 (see also Fig. 11).
In the text

	Fig. 10. As Fig. 9: result obtained with the MIST isochrones for the whole sample. The best fit is the isochrone of 3.8 Myr with an AV of 1.5 (see also Fig. 11), with a relatively large range in age.
In the text

Fig. 11. The best fit MIST and PARSEC model isochrones to the full population, the main sequence only, and pre-main sequence only population, respectively. The absolute magnitude MG is based on the EDR3 parallax. It is clear that the isochrones that was fitted to the pre-main sequence only miss the fit due to the absence of the constraint from the MS stars, resulting in inconsistent extinction, regardless of models. Of the 4 remaining isochrones, the range in the age is between 3.8 Myr and 4.8 Myr and the AV that was determined is around 1.5–1.7 for all models, consistent with the literature values.

In the text

	Fig. 12. Distribution of solutions based on the bootstrap method. The location of the circle represents the age and extinction of the best fit isochrone, and the size of the cirle is proportional to the occurrence rate based on the 200 trials.
In the text

	Fig. 13. Location of HD 153919/4U 1700-37 and the individual members of NGC 6231 now (blue) and 2.2 Myr ago (green). The orange arrows indicate the movement of HD 153919 and NGC 6231 in 0.5 Myr, respectively.
In the text

	Fig. 14. The position of HD 153919 and that of the individual members of NGC 6231 now (blue) and 2.2 Myr ago (green) in Galactic Cartesian coordinates. The uncertainty in the distance of HD153919 is represented by the magenta lines. The orange arrows show the movement of HD 153919 and NGC 6231 in 0.5 Myr.
In the text