The Mass Function of GX 339–4 from Spectroscopic Observations of Its Donor Star^*

We obtained 16 VLT/X-shooter spectra of GX 339−4 in service mode while the system was in quiescence: 1 pilot observation on 2016 May 22 (run code 097.D-0915) and subsequently 15 observations between 2016 August 10 and September 8 granted as director's discretionary time (run code 297.D-5048; see Table 1 for a log of the observations). Every spectrum was obtained in one observing block (OB) including the target acquisition and one ABBA nod pattern, using the XSHOOTER_slt_obs_AutoNodOnSlit template and a nod throw of 5''. Individual exposure times are 250, 200, and 275 s for the UVB, VIS, and NIR arms, respectively. This yields total integration times of 1000, 800, and 1100 s for the respective arms for each of the 16 observations. We use the 10 slit in the UVB arm, the 09 slit in the VIS arm, and the 0.9JH" slit (with the K-band blocking filter) in the NIR arm. This setup was chosen to maximize the signal-to-noise ratio (S/N) in the J- and H-bands. The slit width-limited resolution $R=\lambda /{\rm{\Delta }}\lambda$ is 4350 in the UVB arm, 7410 in the VIS arm, and 5410 in the NIR arm. On some observing nights, more than one OB was executed. In these cases the target was reacquired at the beginning of each OB. A telluric standard star was observed for every science spectrum. Sky conditions were photometric on 2016 August 14 and clear on all other nights.

For all observations we use a sky position angle of $119^\circ$ so that both the counterpart to GX 339−4 and a nearby star (star B in Shahbaz et al. 2001) fall in the slits of the three arms. The separation between GX 339−4 and star B is only $1\buildrel{\prime\prime}\over{.} 7$, so we require a seeing of $\lt 0\buildrel{\prime\prime}\over{.} 9$ (as measured in the V-band) to separate the two spectra. This condition was met in most observations. Because in the NIR arm the seeing is smaller than the slit width in all observations, the resolution is set by the seeing rather than by the slit. This means that the effective resolution varies between the spectra and is somewhat higher than the slit-limited resolution of R = 5410, or ∼56 km s⁻¹ at 1.65 μm. To quantify this, we measure the full width at half maximum (FWHM) of the spectra by fitting two Gaussians to the spatial profile of star B and GX 339−4 at several wavelengths (see Table 1). The spectra taken under the best seeing conditions have a FWHM in the H-band of $\sim 0\buildrel{\prime\prime}\over{.} 65$, which corresponds to a resolution of R = 7490 or ∼40 km s⁻¹.

As expected, there is hardly any signal in the UVB arm. The S/N in the VIS arm is rather low ($\lesssim 5$) but the Hα line is clearly detected. The NIR arm has the highest continuum S/N in all observations, varying from $\sim 5\mbox{--}15$ depending on the seeing. We analyze the data from the VIS and NIR arms only.

We process the data of GX 339−4 and the telluric standards using the X-shooter pipeline in the Reflex environment (Freudling et al. 2013). This produces a flat-fielded, sky-subtracted, rectified, wavelength- and flux-calibrated 2D image and a 1D extracted spectrum. We also use the Starlink package Figaro to optimally extract the GX 339−4 spectra with the tasks profile and optextract (Horne 1986), using a non-symmetric extraction region to avoid as much contamination as possible from star B. Comparing the spectra extracted in these different ways, we find that the optimal extraction significantly reduces the errors for the NIR arm spectra, while it does not improve the VIS arm spectra that are dominated by emission lines. Therefore, we use the optimally extracted spectra in the NIR arm and the spectra provided by the pipeline for the VIS arm. For the telluric standards we take the normally extracted spectra provided by the pipeline.

The size of the extraction aperture for the NIR spectra is determined by the proximity of star B, not by the width of the spatial profile. This causes a loss of NIR flux that is most significant in the spectra taken under the worst seeing conditions. In combination with the fact that most spectra were taken under clear, not photometric conditions, this means that the absolute flux calibration of our spectra is not accurate.

For the telluric correction of the NIR arm data, we use the telluric calibration package Molecfit (Kausch et al. 2015; Smette et al. 2015). First, we fit an atmospheric model to the telluric standard stars, to determine the concentrations of H₂O, CO₂, CH₄, and O₂. We then fit an atmospheric model to the GX 339−4 spectra, keeping the concentrations fixed, and apply the telluric correction.

To measure the RVs and rotational broadening of the GX 339−4 spectra we compare them with template spectra. We use both spectra of Galactic stars observed with X-shooter, and synthetic spectra from model atmospheres at different temperatures, as templates.

We retrieved the raw data of several K-stars and one late G-type star from the ESO archive (see Table 2). These stars were observed with X-shooter as part of the X-shooter Spectral Library (XSL) project, with the 09 NIR arm slit, as used for our GX 339−4 data, but without the K-band blocking filter (Chen et al. 2014, program ID 085.B-0751). All were observed with a seeing of $\gt 0\buildrel{\prime\prime}\over{.} 9$ in the NIR, so their spectral resolution is determined by the slit width. The XSL contains only one K-type sub-giant; to obtain templates spanning a range of spectral types we also include several giant stars. We reduce the NIR arm data the same way we reduced that for the telluric standard stars, using the X-shooter pipeline in Reflex to extract the spectra and Molecfit to correct for telluric absorption.

In addition, we use three model atmospheres with $\mathrm{log}g=2.5$, solar metallicity, and effective temperatures of 3938, 4500, and 5000 K. These were produced using the state-of-the-art realistic three-dimensional radiative hydrodynamical (RHD) simulations of stellar convection carried out with the Stagger-code (Nordlund et al. 2009). These simulations cover a large part of the HR diagram, including the evolutionary phases from the main sequence over the turnoff up to the red-giant branch for low-mass stars (Magic et al. 2013). These simulations have been used to compute synthetic spectra in the range of X-shooter with the multidimensional pure-LTE radiative transfer code Optim3D (Chiavassa et al. 2009). The code takes into account the Doppler shifts occurring due to convective motions and solves monochromatically for the emerging intensity, including extensive atomic and molecular continuum and line opacity data from UV to far-IR. We smooth the high-resolution RHD simulations both to the slit width-limited X-shooter resolution (56 km s⁻¹) and to the highest resolution present among our GX 339−4 observations (40 km s⁻¹) and rebin the two sets of synthetic spectra to the same velocity scale as the data.

We use Tom Marsh's program Molly⁷ to further analyze the data. First, we apply hfix to convert to heliocentric velocities and we mask noisy pixels (mainly due to residuals from strong telluric emission lines). For the NIR data, we use vbin to rebin all spectra to a common velocity spacing of 15 km s⁻¹ per pixel, the average velocity dispersion of X-shooter in the J- and H-bands. Absorption lines corresponding to Mg i ($1.71\,\mu {\rm{m}}$) and Al i (at $1.31\,\mu {\rm{m}}$) are visible in the individual spectra, indicating that the donor star is a K-type star (e.g., Meyer et al. 1998; Wallace et al. 2000, see Section 3.2). To prepare the spectra for cross-correlation, we first scale them by dividing through a constant and then subtract a spline fit. K-star spectra have many absorption lines, which makes defining the continuum level difficult. We determine the normalization factor by fitting a constant to a relatively clean region of the spectrum close to the strong Mg i line at 1.71 μm. For the spline fit, we mask the regions in the J-band and between the J- and H-bands that are severely affected by telluric absorption, as well as the H and He emission lines that are present in the GX 339−4 spectra (see Figure 1).

Zoom In Zoom Out Reset image size

We prepare the templates the same way we prepared the GX 339−4 spectra. In addition, we cross-correlate the spectra of the templates stars observed with X-shooter with the RHD simulations to obtain accurate radial velocities, and shift the template spectra to remove their intrinsic radial velocity. We then cross-correlate the 16 individual GX 339−4 spectra with the different templates using the task xcor, excluding regions with high noise levels and emission lines. All three RHD simulations, at both resolutions, and the XSL stars, yield consistent RVs. The highest values for the cross-correlation function are found for the T = 3938 K model smoothed to 56.6 km s⁻¹, hence the results in this paper are from the cross-correlation with that template. For five GX 339−4 spectra with lower S/N we find that only considering the regions around the strongest absorption lines in the cross-correlation significantly improves the signal. The highest-S/N spectra, and low-S/N spectra with many masked pixels in these regions, do not benefit from this approach—for those we use the RV found by cross-correlating the full selected spectral range.

xcor gives only a purely statistical error that underestimates the real uncertainty in the RV measurement. To better constrain the uncertainty we use the boot command in Molly to produce 1000 bootstrapped copies of every spectrum, which we then cross-correlate with the model spectrum. We fit a Gaussian to the resulting distribution of RVs, and adopt the mean and standard deviation of that Gaussian as reliable estimates of the RV and uncertainty. The results are listed in Table 1; the 1–σ error includes the uncertainty in the wavelength calibration of X-shooter (2 km s⁻¹).

To determine the orbital period we fit sinusoids of the form $\gamma +{K}_{2}\sin (2\pi (t-{T}_{0})/P)$ for a range of periods and calculate ${\chi }^{2}$ (see Figure 2). We find that the minimum ${\chi }^{2}$ is 25.6 (with 13 degrees of freedom) for a period of $P=$ 1.7587 ± 0.0005 days, where the 1σ uncertainty is set by the change in P for which ${\rm{\Delta }}{\chi }^{2}=1$. Our spectroscopic orbital period confirms the orbital period reported by Hynes et al. (2003) and Levine & Corbet (2006). The resulting RV curve is plotted in Figure 3. The best-fit parameters are γ = 26 ± 2 km s⁻¹, measured with respect to the heliocentric velocity, K₂ = 219 ± 3 km s⁻¹, and the time of inferior conjunction of the donor star ${T}_{0}=$ MJD 57529.397 ± 0.003.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, we evaluate K₂ using the relation between this orbital parameter and the FWHM of the Hα line (Casares 2015). The Hα line in our GX 339−4 spectra is double-peaked. Following the approach of Casares (2015), we fit a single Gaussian profile to the Hα emission line in each VIS spectrum with mgfit. We find an average FWHM of 880 ± 20 km s⁻¹(corrected for the instrumental resolution; see Figure 4). This corresponds to ${K}_{2}=206\pm 12$ km s⁻¹ according to the Casares (2015) relation.

Zoom In Zoom Out Reset image size

To obtain a higher-S/N spectrum for determining both the projected rotational velocity and spectral type of the donor star, we normalize the individual NIR spectra by dividing them through a spline fit, correct for the RVs of the individual spectra, and average them, weighted by their mean (S/N)². In the resulting spectrum, several absorption lines are detected, mainly due to Al i and Mg i (see Figure 5). The lines are weaker than those in the model spectra, indicating that the accretion disk still contributes significantly to the detected NIR continuum flux, and is broadened because of the rotation of the tidally locked star.

Zoom In Zoom Out Reset image size

To determine the projected rotational velocity $v\sin i$ of the donor star, and the contribution of the accretion disk to the NIR emission, we use the optimal subtraction method (Marsh et al. 1994). In Molly, we broaden all templates by a range of velocities (0—100 km s⁻¹ in steps of 1 km s⁻¹) with rbroad and then use optsub, with a linear limb darkening coefficient of 0.75, to subtract the broadened models from the average GX 339−4 spectrum. For the optimal subtraction we use small spectral regions containing the five strongest absorption lines in the GX 339−4 spectrum. optsub subtracts a constant times the template from the GX 339−4 spectrum, adjusting the constant to minimize the residual scatter between the spectra. The scatter is measured by computing the ${\chi }^{2}$ between the residual spectrum and a smoothed version of itself, produced by convolving the subtracted spectrum with a Gaussian with a FWHM of 25 pixels.

Apart from the two worst-matched templates, all ${\chi }^{2}$ distributions yield consistent $v\sin i$ values, with only small differences between the results found using the higher- and lower-resolution RHD simulations (see Figure 6). This confirms that we are resolving the stellar absorption lines with our instrumental setup under the prevailing seeing conditions. We find that using a different limb darkening factor of 0.5 influences the $v\sin i$ measurement by $\lt 5$ km s⁻¹. The best-fitting templates are the T = 3938 K RHD simulation, a K2III star and the K1IV star; the K2III star HD 175545 has the lowest minimal ${\chi }^{2}$ of 278 (with 361 degrees of freedom), for an optimum factor of 0.47 ± 0.04. This template, and the residuals of the optimal subtraction, are shown in Figure 5. To obtain a reliable estimate of the uncertainty we follow the bootstrapping approach by Steeghs & Jonker (2007): we produce 1000 bootstrapped copies of the average GX 339−4 spectrum and repeat the optimal subtraction analysis with the 2 best-matched templates (HD 175545 and HD 165438) for each of them. For every bootstrapped copy we find the value of $v\sin i$ that results in the minimum ${\chi }^{2}$. We fit a Gaussian to the resulting distributions of ${\chi }^{2}$ values and adopt the mean and standard deviation as the best estimate of $v\sin i$ and its uncertainty. For both templates we find that $v\sin i$ = 64 ± 8 km s⁻¹, and the donor star contributes $\sim 45 \% \mbox{--}50 \%$ of the light in the NIR.

Zoom In Zoom Out Reset image size

Our measurement of K₂ is significantly lower than the lower limit of 317 km s⁻¹ reported by Hynes et al. (2003) based on the Bowen blend emission lines. One explanation for this discrepancy could be that the absorption lines in our NIR spectra are predominantly formed very close to the tip of the Roche lobe, and our K₂ measurement simply requires a larger K-correction than the Bowen blend measurement. Following the procedure for K-correction outlined by Muñoz-Darias et al. (2008), taking the maximum acceptable value for f with $q\approx 0.18$ based on our $v\sin i$ measurement (see Section 4.2), we find a maximum correction factor of ∼1.7. This could in principle reconcile our value for K₂ with the one found by Hynes et al. (2003). This scenario would imply significant X-ray irradiation of the donor star during quiescence and that the absorption lines are stronger when we view the irradiated side of the donor star. Our spectra cover most orbital phases; we find no such variation in the presence or strength of the absorption lines. In contrast, the Bowen blend emission lines used by Hynes et al. (2003) are only present in part of the orbit. The total NIR flux does vary between our spectra, but not as a function of phase; it is correlated to the observing conditions. This is because the size of the extraction region is limited by the close proximity of the second star in the slit rather than by the width of the profile, which causes more flux to be lost in spectra taken under worse seeing conditions. These results support that our measured RVs reflect the motion of the center of mass of the donor star and we do not underestimate K₂.

In addition, the K₂ implied by the FWHM of the Hα emission line is 206 ± 12 km s⁻¹, consistent with the value obtained from our dynamical study. As noted by Casares (2015), a caveat when using this relation is that the width of the line rises steadily after an outburst and it takes several years of quiescence to "settle down." Because of the frequent outbursts of GX 339−4 it is not clear if this system ever reaches that equilibrium state. For comparison, Rahoui et al. (2014) obtained three spectra of GX 339−4 during its 2010 outburst, two while it was rising in the hard state and one after it transitioned to the soft state, and report values for the FWHM of the Hα line of 550 ± 30, 704 ± 53, and 421 ± 44 at those epochs, respectively. Our FWHM of 880 ± 20 is significantly higher, in line with expectations. The last major outburst of GX 339−4 took place in 2014 December, 1.5 years before our observations. In V404 Cyg, for which Casares (2015) shows the evolution of the Hα FWHM in the 20 years after an outburst, the width of the line approaches its equilibrium state after ∼500 days (their Figure 2). This suggests that our K₂ measurement from this relation should be close to the real value, which further strengthens the case for a lower K₂ than that found through the Bowen blend method.

This large discrepancy between the RV semi-amplitude measured from the Bowen blend lines and our value shows that in the case of GX 339−4, the assumption that the narrow emission lines originate in the irradiated side of the donor star is likely wrong. Rather, the sharp Bowen blend lines must predominantly be emitted in a different location in the system with a higher Keplerian velocity or be a combination of the Keplerian motion of the star and gas flow motions. This implies that using the RV curve of the narrow Bowen blend lines may not always be a good alternative to obtaining the RV curve of the donor star by the motion of photospheric absorption lines if one wants to measure K₂ in BH XRBs.

Our K₂ measurement corresponds to a mass function

$$\begin{eqnarray}&&f(M)=\displaystyle \frac{{K}_{2}^{3}{P}_{\mathrm{orb}}}{2\pi G}=\displaystyle \frac{{M}_{x}{\sin }^{3}i}{{(1+q)}^{2}}=1.91\pm 0.08\,{M}_{\odot },\end{eqnarray} \tag{ 1 }$$

where M_x is the mass of the compact object and q is the mass ratio between the donor star and the compact object, ${M}_{\mathrm{donor}}/{M}_{x}$. This mass function is a factor ∼3 lower than that derived from the Bowen emission line measurements for GX 339−4.

The mass ratio q in a Roche lobe overflow system is related to the projected rotational velocity and K₂ as

$$\begin{eqnarray}&&\displaystyle \frac{v\sin i}{{K}_{2}}=0.462{q}^{1/3}{(1+q)}^{2/3},\end{eqnarray} \tag{ 2 }$$

given that the rotation of the mass donor star will be tidally locked with the orbital period (e.g., Gies & Bolton 1986). We find that $v\sin i=$ 64 ± 8 km s⁻¹, which, combined with our measurement of K₂, implies that q = 0.18 ± 0.05.

The inclination of GX 339−4 has been suspected to be low rather than high based on the small separation of the double-peaked emission lines (Wu et al. 2001).

Modeling of the reflection component of the X-ray spectra of GX 339−4 has yielded inclinations ranging from 30° to 60° (e.g., Fürst et al. 2015; García et al. 2015; Basak & Zdziarski 2016; Parker et al. 2016). An inclination in this range would lead to a BH mass between $\sim 4\mbox{--}20$ M_⊙ (see Figure 7). However, this method gives the inclination of the inner accretion disk, which is not necessarily the same as the binary inclination if the BH spin angular momentum axis is not aligned with the binary angular momentum vector.

Zoom In Zoom Out Reset image size

An upper limit to the binary inclination is provided by the fact that GX 339−4 does not show eclipses in its X-ray light curve; with q = 0.18 ± 0.05, the donor star subtends an angle of $\sim 25^\circ$, as seen from the compact object (Paczyński 1971), setting an upper limit to the inclination of $i\,\lt$ 78°. Combined with our q value, at this inclination our mass function yields a BH mass of 2.9 ± 0.3 M_⊙, and a 95% confidence lower limit to the mass of the compact object of ${M}_{x}\,\gt$ $2.3\,{M}_{\odot }$.

Muñoz-Darias et al. (2008) found a maximum mass for a sub-giant donor star of 1.1 M_⊙. Combined with our mass ratio, this implies a maximum mass for the BH of ${M}_{x}=6.1\,\pm 1.7$ M_⊙, or a 95% confidence upper limit of $9.5\,{M}_{\odot }$. This also sets a lower limit to the binary inclination of 37°.

A reliable measurement of the binary inclination of GX 339−4, for example, through the detection of ellipsoidal modulations, is required to determine the mass of the accretor in this system.

From our optimal subtraction analysis we conclude that the donor is most likely a K1 or K2 star, confirming the results of Hynes et al. (2004) and Muñoz-Darias et al. (2008). The S/N of our spectrum is not high enough to significantly detect many of the fainter absorption lines that are present in such stars, but we clearly detect the strongest lines due to Mg i and Al i(see Figure 5). However, the strongest Si i line in the template spectra (at 1.59 μm) is not detected in our GX 339−4 spectrum, which may indicate that the donor has an anomalous Si abundance. Higher-S/N spectra are necessary to investigate this further. The star contributes $\sim 45 \% \mbox{--}50 \%$ of the NIR emission.

The distance to GX 339−4 is not accurately known; Hynes et al. (2004) found a lower limit to the distance of 6 kpc, but point out that their observations are also consistent with a much larger distance of ∼15 kpc. The flux calibration of our spectra is not absolute, as they were largely taken under non-photometric conditions and because the size of the extraction region is limited by the proximity of star B. However, we can derive a lower limit to the distance of GX 339−4 using the spectrum with the highest measured flux. This is the spectrum taken on 2016 May 22, which has an average flux density of $1.2\times {10}^{-17}$ erg cm⁻² s⁻¹ Å⁻¹ in the H-band—we use the H-band because it suffers the least from extinction. Taking 50% as the stellar contribution in line with our optimal subtraction results, this translates into an H-band magnitude of the donor star of 18.3.

For a lower limit on the distance we take the minimum mass models for the donor star from Muñoz-Darias et al. (2008), the properties of which translate into an absolute H-band magnitude of ∼2.8. For the extinction toward GX 339−4 we follow Muñoz-Darias et al. (2008) and adopt $E(B-V)\approx 1.2$ (Hynes et al. 2004), which corresponds to ${A}_{H}\approx 0.7$. Using these values we find a distance to GX 339−4 of ∼9 kpc. The main source of uncertainty is in the flux of the donor star. We conservatively estimate the uncertainty on the stellar flux as 50%. This translates into a 25% uncertainty on the distance. We thus find a lower limit on the distance of ∼5 kpc, although given the fact that the X-ray activity of GX 339−4 favors a donor mass at the higher end of the allowed mass range (Muñoz-Darias et al. 2008), a larger distance seems more likely. These conclusions are in agreement with those of Hynes et al. (2004).

It is interesting to note the difference in K₂ between GX 339−4 and XTE J1550-564, which has a very similar orbital period (1.54 days) and has been suggested as a "twin" system to GX 339−4 (Muñoz-Darias et al. 2008). Its donor star has been identified as a K3 star (Orosz et al. 2011). However, XTE J1550-564 in quiescence shows an Hα line with a FWHM of 1506 ± 151 km s⁻¹ (Casares 2015), and it has an RV semi-amplitude of 363.14 ± 5.97 km s⁻¹ (Orosz et al. 2011). The difference between the two systems may be due to inclination: XTE J1550-564 has a relatively high inclination of $\sim 75^\circ$, and an inclination of $\sim 35^\circ$ for GX 339−4 would explain the difference in K₂. Alternatively, GX 339−4 may contain a less massive BH than XTE J1550-564, which contains a BH with a mass in the range of 7.8–15.6 M_⊙. The lower mass ratio in XTE J1550-564 ($q\approx 0.03$) would suggest a more massive BH if the donor star masses are similar in both systems.

We obtained 16 VLT/X-shooter observations of GX 339−4 and detected absorption lines from the donor star in the NIR spectrum, allowing us to measure the RV semi-amplitude and projected rotational velocity of the donor star for the first time. We measure K₂ = 219 ± 3 km s⁻¹ and $v\sin i$ = 64 ± 8 km s⁻¹. This implies a mass ratio q = 0.18 ± 0.05 and a mass function f(M) = 1.91 ± 0.08 ${M}_{\odot }$. The value we find for K₂ is significantly lower than the value obtained using the Bowen blend method, showing that a mass function based on RV measurements of these emission lines can be incorrect.

The donor is a K1-2 type star and we estimate that it contributes $\sim 45 \% \mbox{--}50 \%$ of the light in the J- and H-bands. Adopting the minimum mass model for a stripped donor star from Muñoz-Darias et al. (2008), we obtain a lower limit to the distance to GX 339−4 of ∼5 kpc, in agreement with Hynes et al. (2004).

Without the binary inclination we can only set limits to the accretor mass; we find $2.3\,{M}_{\odot }$ $\leqslant \,{M}_{{\rm{X}}}\,\leqslant$ $9.5\,{M}_{\odot }$, which means that a (massive) neutron star accretor cannot be excluded. Although a low system inclination would still allow for a BH more massive than 5 M_⊙, the mass function of GX 339−4 is much lower than has been assumed to date and it may in fact be the first BH to fall in the "mass-gap" of 2–5 M_⊙ (Özel et al. 2010; Farr et al. 2011).