DISTRIBUTION OF HIGH-MASS X-RAY BINARIES IN THE MILKY WAY

Our fitting procedure is based on the Levenberg–Marquardt least-squares algorithm implemented in Python (leastsq routine from the scipy package). For each HMXB, we build the SED in optical and NIR from a maximum of 8 (U, B, V, R, I, J, H, Ks) and a minimum of 4 mag points. This SED is then fit (see Figure 1 for an example) by a blackbody model, with the aim of using a uniform fitting method for the entire sample of HMXBs, given by the relation

$$\begin{equation} \lambda F_{\lambda } =\frac{2\pi hc^2}{\lambda ^4}\times 10^{-0.4A_{\lambda }}\frac{\left({R/D}\right)^2}{\exp {\left(\frac{hc}{\lambda k_{B}{T}}\right)}-1}, \end{equation} \tag{ 1 }$$

where λ is the wavelength in μm, F_λ is the flux density in W m⁻² μm⁻¹, h is the Planck constant, c is the speed of light, A_λ is the extinction at the wavelength λ, R/D is the stellar radius over distance ratio, k_B is the Boltzmann constant, and T is the temperature of the star. When the radius and the temperature of the companion star are available in the literature, we use these values (given in Table 2). For the other sources, we derive the radius and the temperature of the companion star, which dominates the optical and NIR flux, from the spectral type and the luminosity class, by using tables of Vacca et al. (1996), Panagia (1973), Martins & Plez (2006), and Searle et al. (2008), hereafter noted PVMS. Since only four sources (including GX 301−2) with a highly evolved mass star have been detected in the Milky Way (Mason et al. 2012), we can reasonably assume that for all the other systems studied here, the radius of the companion star is sensibly close to the one expected, depending on the luminosity class. When the spectral type is poorly known, we derive the mean distance value by considering the different possible spectral types and allocate adequate error on the distance of these systems. Finally, if no data are given in Vacca et al. (1996), Panagia (1973), and Searle et al. (2008) about a certain spectral type (which is the case for luminosity classes II and IV), we computed intermediate R and T between luminosity classes I and III for class II and between classes III and V for class IV, respectively. This affects only three sources in our sample. For stars with peculiarities (N iii and He ii in emission mentioned by an "f" in the spectral type, broad absorption mentioned by an "n," unspecified peculiarity mentioned by a "p") we assume their radius and temperature to have the same values as "normal" stars. For Be stars, we take R and T values of normal B stars, and the circumstellar disk emission is taken into account through a possible uncertainty due to the expected infrared excess of this family of sources (see Section 2.2). Finally, we describe our estimate of distance uncertainties due to errors on R and T of the companion star in Section 2.2.

Zoom In Zoom Out Reset image size

Moreover, two parameters are left free: the extinction in V band A_V and the ratio R/D, whereas the extinction A_λ is derived at each wavelength from Cardelli et al. (1989) assuming R_V = 3.1 (corresponding to the mean value derived for a diffuse interstellar medium toward the Galactic plane; Fitzpatrick 1999). The influence of potential change of the extinction law was explored. The results, presented in Section 2.2, do not show any substantial variation in distance determination due to extinction change.

Knowing the radius R of the companion star, we then calculate the distance D in kpc.

The least-squares function given by the formula χ² = ∑_i[(X_{i, obs} − X_{i, model})/(σ_i)]² (with X_{i, obs} the observed flux value for the filter i, X_{i, model} the predicted flux in the ith filter derived from the blackbody model, and σ_i the flux error in the same filter) is then minimized by the Levenberg–Marquardt algorithm.

To check our results, we estimated the distance of each HMXB using the expression

$$\begin{equation} D_{\rm {pc}} = 10^{-0.4(m_{\rm {V}}-M_{\rm {V}}-A_{\rm {V}}+5)}, \end{equation} \tag{ 2 }$$

with m_V the relative magnitude and M_V the absolute magnitude while A_V is given by the fit. This formula has the advantage of not depending on the companion star radius. However, it does depend on the absolute magnitude M_V, given in Martins & Plez (2006) for O-type stars and in Morton & Adams (1968) and Panagia (1973) for B stars. The results obtained with this second method are consistent with those using the first approach.

We give results in Table 1 and compare these results with previously published results in Table 3. We point out that median discrepancy in distance is ∼17% and A_V is often very similar (median discrepancy of less than 7%).

The magnitude uncertainties are retrieved from the literature. For sources for which no error is given, we use a systematic error of 0.1 mag, conservative of typical errors. The flux uncertainties are then derived from these magnitude errors. Nonetheless, we assume in the following that the spectral type given in the literature is the real spectral type of the companion star. We carried out simulations on a wide range of spectral types, which enabled us to constrain the uncertainty that appears to be less important for supergiant stars than for zero-age main-sequence stars.

Uncertainties on the radius and temperature of the companion star could have a severe influence on distance and A_V error. We generate SEDs of various sources using a large range of temperatures, which clearly show that the distance error is less sensible to the temperature error than to the radius one. That can be easily understood given the fact that distance is directly derived using the R/D ratio. We distinguished different cases. (1) For sources for which spectral type is accurately determined, we select R and T values mostly in Vacca et al. (1996) for consistency. If no data are given in this reference, we used R and T values given in Panagia (1973), Martins & Plez (2006), or Searle et al. (2008). With the aim of deriving accurate distance and A_V errors due to radius and temperature uncertainties, we fitted the SED of each system using the R and T values given in PVMS, and we compute the dispersion in distance and A_V values obtained using these four references. These computations lead to a mean difference in distance value of 10%, as well as for the A_V value. Then, for systems for which spectral type is well known, we consider a systematic error of 10% of the derived distance and A_V values due to radius and temperature uncertainties. (2) For sources for which values of R and T are available in the literature (see Table 2), we use these values in our fitting procedure to derive distance and extinction in V band. When maximum and minimum radius and temperature values are available, we derive distances and dispersion on distance using these values. We also derive the distance considering R and T values given in PVMS. Then, we consider as the final error (due to radius and temperature uncertainties) the dispersion on all the distance values obtained with R and T values given in the four references and with R and T values given in the paper dedicated to the source (referenced in Table 2). The same approach was used to derive an A_V error due to radius and temperature uncertainties. (3) When the spectral type is not well constrained and if no R and T value is available in the literature, we derived distances (and then a mean distance) using the different possible spectral types of the companion star and then the different R and T values given in PVMS. We finally compute the dispersion on all the distances derived and considered it as the error due to radius and temperature uncertainties. We derive an A_V error due to radius and temperature uncertainties for these sources, using the same method. (4) For luminosity classes II and IV, since no data are given in PVMS, we used intermediate R and T values between classes I and III for class II and between classes III and V for luminosity class IV. Errors were computed using the dispersion on distance obtained with these R and T values and distances obtained with R and T values of lower and higher luminosity classes.

Degeneracy between several parameter values (based on the fitting procedure) needs to be taken into account. Indeed, solely relying on a single best fit does not capture the full phenomenology associated with SED fitting because D and A_V are degenerate in this approach. In order to produce the best set of fits and to determine the dispersion on distance and extinction, we carried out 500 Monte Carlo simulations for each observed source, by varying the photometry within the uncertainties. Hence, we generate a random number from a normal distribution (assuming the photometric errors to be Gaussian), defined by the error bars for each photometric point, so that we build 500 new SEDs derived from the original one. These 500 new SEDs are subject to the same χ² statistic computation as the one described above. Then, we have an entire set of best fits of parameters (D, A_V) and are able to plot the distribution in the parameter space, showing the distribution of properties derived from these Monte Carlo simulations, and especially showing the dispersion on distance and extinction for each source (see Figure 2). This dispersion value is taken as the error of the fitting procedure, and the median value of dispersion on distance determination for all considered HMXBs is then 0.75 kpc.

Zoom In Zoom Out Reset image size

There are other sources of uncertainties, particularly the infrared excess of Be stars due to their circumstellar envelope generating free–free radiation. According to Dougherty et al. (1994), this excess should not exceed a mean 0.1 mag in J band, 0.15 mag in H band, and 0.25 mag in K band. However, this value corresponds to absolute magnitude, and therefore the excess can be smaller in apparent magnitude for sources located far away and higher for close ones. To take this effect into account, an estimate of the distance and extinction is needed. Since these two values are derived from the fitting procedure, we are only able to consider the distance and absorption values obtained without taking this IR excess into account. This approach is finally equivalent to adding a conservative error of 0.1, 0.15, and 0.25 mag to the apparent magnitude in the J, H, and Ks bands, respectively; this method presents the advantage of treating all sources in a uniform way. Based on these results, presented in Figure 3, we assume that this uncertainty will not affect the distances derived in this paper.

Zoom In Zoom Out Reset image size

Furthermore, to test the potential additional error due to a different extinction law, we use R_V values given by Geminale & Popowski (2004). We affected to each HMXB an R_V value derived from the closest star present in their catalog, and we fitted the SED using the Cardelli law with the new R_V value. We find a median difference in distance determination of 0.07 kpc. Since deriving an accurate R_V value for each system appears to be observationally biased (it depends on the closest star available in their catalog), we decided not to take into account this additional uncertainty.

Finally, Table 1 presents the 46 sources and their fundamental parameters (computed in this work or taken in the literature). Figure 10 represents most of the studied HMXBs with the uncertainties on their location computed taking into account all the errors described above: the error on distance due to uncertainties on the companion star radius and temperature and the error coming from the fitting procedure. We derived in the same way the final error on A_V.

To perform this study, we need to consider a number of different issues. As underlined above, the expected offset depends on the age of the X-ray sources; hence, to highlight this offset, we must split the sample of HMXBs depending on the age of the sources. Two different samples are then created: one containing four supergiant stars (luminosity class I or II according to Charles & Coe 2006) and a second one containing nine Be stars (luminosity class III or V according to Charles & Coe 2006). We explain the way these samples were created in Section 4.2.

We calculate the distance from each HMXB to the closest actual spiral arm (given in Russeil 2003) using

$$\begin{equation} \rm {distance}=\sqrt{(\it X_{\rm {source}}-X_{\rm {arm}})^2+(Y_{\rm {source}}-Y_{\rm {arm}})^2}, \end{equation} \tag{ 5 }$$

where (X_source, Y_source) are the coordinates of the source and (X_arm, Y_arm) are the coordinates of the closest point on the arm.

We follow exactly the same procedure to calculate the distance from each source to the closest expected position of sources formed 20, 40, 60, 80, and 100 Myr ago, and we determine the mean value of the offsets (taking into account all the sources of the two samples). The sources are expected to be closer to one of the expected positions computed above than to the current spiral arms observed by Russeil (2003). The method is described in Figure 11 (for instance, in Figures 10 and 11, a source of 20 Myr should be located closer to the green arm representing the expected position of a 20 Myr old HMXB than to the current spiral wave position). Results are given in Figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We determine the error bars using the propagation of uncertainty formula. The 1σ error associated with the distance between each source i and the closest point on the arm is then given by the following equation:

$$\begin{eqnarray} \sigma _{i}^2 = & \frac{(X_{\rm {source}}-X_{\rm {arm}})^2\left[\sigma _{X_{\rm {source}}}^2+\sigma _{X_{\rm {arm}}}^2\right]}{(X_{\rm {source}}-X_{\rm {arm}})^2+(Y_{\rm {source}}-Y_{\rm {arm}})^2} \nonumber\\ \ms & +\frac{(Y_{\rm {source}}-Y_{\rm {arm}})^2\left[\sigma _{Y_{\rm {source}}}^2+\sigma _{Y_{\rm {arm}}}^2\right]}{(X_{\rm {source}}-X_{\rm {arm}})^2+(Y_{\rm {source}}-Y_{\rm {arm}})^2}, \end{eqnarray} \tag{ 6 }$$

where ($\sigma _{X_{\rm {source}}}$, $\sigma _{Y_{\rm {source}}}$) are the errors on source position detailed in Section 2.2 and ($\sigma _{X_{\rm {arm}}}$, $\sigma _{Y_{\rm {arm}}}$) are the errors on arm position taken to be equal to zero. Finally, the uncertainty associated with the mean value is calculated as follows:

$$\begin{equation} \sigma =\frac{1}{N}\sqrt{\sum _{i}^{N}\sigma _{i}^2}, \end{equation} \tag{ 7 }$$

where N is the number of sources and σ_i is the error associated with the ith source.

Taking into account all the sources of the sample (Figure 12, lower panel), we cannot observe any significant variation of the offset with time before 60 Myr. This result was expected because the minimum peak of distance between each HMXB and the closest expected position should be clearly different whether we consider Be or supergiant stars. Also, by taking all spectral types into account, we tend to lose a part of the information except the mean age upper limit of 60 Myr highlighted by the plot. For both supergiant and Be samples, we observe a more significant increase of the offset after 60 Myr (Figure 12, upper and middle panels). This increase could be a signature of the expected offset between the HMXB positions and the current spiral arms. However, we must be cautious about this result especially because of the small number of HMXBs in the two samples.

Even if the locations of the sources are accurately determined, several reasons may affect the HMXB density and prevent the offset detection as underlined by Lutovinov et al. (2005): complex motion of density wave and stars from their birthdate to the X-ray phase, presence of previously undetected parts of the Galactic spiral arms, observational selection effect, etc. Finally, a larger sample of supergiant type HMXBs is needed to confirm the offset detection more confidently.

Moreover, the time interval during which HMXBs appear should translate the mass range of both stars of binary systems (see Dean et al. 2005). The results presented in Figure 12 only enable us to state that this time interval is lower than 60 Myr on average for all stars.

Following the method described above, one can compute a distance between each source and the different theoretical arms that correspond to the current predicted position of a source sample born 20, 40, 60, 80, or 100 Myr ago. In the previous section we studied a sample of Be and supergiant HMXBs, and this could also be applied to each source separately, to constrain the age and the potential migration of the system due to a supernova kick.

We choose not to take into account here the arm width that could be translated as an uncertainty on the position of the expected position of sources. Indeed, by taking this position dispersion into consideration, 1σ errors considerably increase and prevent any conclusion. Then, we assume all sources to be formed at the central position of the arm.

We expect the distance from the source to the expected position to decrease until the expected position corresponds to the age of the source, and then to increase afterward (see Figure 11). For some sources, the offset change is more complicated and can be explained as follows:

• 1.  
  Belonging to one of the four arms of the Milky Way is not well established when the distance between a source and the theoretical expected positions is seen as always increasing with time.
• 2.  
  For sources located close to the corotation radius, current spiral arms and expected positions of HMXBs of different ages are almost superimposed, explaining the quasi-constant distance to theoretical positions observed for some sources. Since offsets between expected positions of sources of different ages are very small, it can be possible to ascribe an incorrect age to a source located close to an expected position that does not correspond to its real age.
• 3.  
  Sources that underwent more significant kicks could be incorrectly associated with a theoretical position of sources of different age.

Then, for this study, we only take into account sources presenting an expected offset evolution during time, to make the conclusions easier.

We present in Figure 13 results for BeHMXBs and sgHMXBs for which the evolution of the offset as a function of time is consistent with the theory. It is possible to determine the rough age of these sources and a lower and upper limit of kick migration distances. The results are summarized in Table 4 and will be discussed in the following paragraphs.

Zoom In Zoom Out Reset image size

If we consider only the sources studied separately here (sources showing the expected evolution of offset as a function of age), we can derive a mean age of ∼45 Myr for sgHMXBs and 51 Myr for BeHMXBs. This result seems consistent with the distinct evolution timescales of these two kinds of systems, but again, we must be cautious about this result because of the small sample of sources used in this calculation (four supergiant and nine Be systems).

Even if the theory is still poorly understood, it is admitted that a small asymmetry in the way a supernova explodes could make a binary system move due to a substantial increase of its velocity in a particular direction. We computed the cluster size (in the same way as in Section 3) only for Be stars that are expected to have experienced a supernova kick event during their evolution, and we found a value of 0.3 kpc, consistent with the migration distances underlined in Bodaghee et al. (2012). This value also enables us to constrain the time lapse between supernova event and HMXB stage; assuming a regular value of kick velocity of 100 km s⁻¹ (for higher kick velocity, fewer systems should remain bound; Hills 1983) and a maximum migration distance due to the kick event of 0.3 kpc, we reach an upper limit for the time lapse between the supernova explosion and the HMXB step of around 3 Myr.

To improve this study, we focus on the sources selected before (four sgHMXBs and nine BeHMXBs; Table 4). The distance between the object and the closest expected position gives a kick value for each source and a mean value depending on the spectral type of the two samples. Results are presented in Table 4. We derive a mean migration distance of 0.11 kpc for BeHMXBs and 0.10 kpc for sgHMXBs. Again, we should be careful about these mean values because of large error bars and small samples. Moreover, it is important to underline that these derived values only represent a lower limit to the kick migration distance since we cannot take into account the migration distance on galactic latitude given by the kick. Then, we only get the projected migration distance on the Galactic plane.