The orbital period, black hole mass and distance to the X-ray transient GRS 1716-249 (=N Oph 93)

We present evidence for a 0.278(8) d (=6.7 h) orbital period in the X-ray transient GRS 1716-249 (=N Oph 93), based on a superhump modulation detected during the 1995 mini-outburst plus ellipsoidal variability in quiescence. With a quiescent magnitude of r=23.19+-0.15 N Oph 93 is too faint to warrant a full dynamical study through dedicated time-resolved spectroscopy. Instead, we apply the FWHM-K2 correlation to the disc Halpha emission line detected in Gran Telescopio Canarias spectra and obtain K2=521+-52 km/s. This leads to a mass function f(M)=4.1+-1.2 Msun, thus indicating the presence of a black hole in this historic X-ray transient. Furthermore, from the depth of the Halpha trough and the quiescent light curve we constrain the binary inclination to i=61+-15 deg, while the detection of superhumps sets an upper limit to the donor to compact star mass ratio q=M2/M1<=0.25. Our de-reddened (r-i) colour is consistent with a ~K6 main sequence star that fills its Roche lobe in a 0.278 d orbit. Using all this information we derive a compact object mass M1=6.4+3.2-2.0 Msun at 68 per cent confidence. We also constrain the distance to GRS 1716-249 to 6.9+-1.1 kpc, placing the binary ~0.8 kpc above the Galactic Plane, in support of a large natal kick.


INTRODUCTION
The Milky Way contains an estimated population of ≈ 10 3 − 10 4 stellar-mass black holes (BH) accreting from a late-type companion star (Yungelson et al. 2006;Kiel & Hurley 2006).As a result of their low average mass-accretion rates these BH binaries are X-ray transients (XRTs) and are routinely discovered at a rate of ∼ 1.7 yr −1 during month-long episodes of increased accretion activity (McClintock & Remillard 2006;Belloni et al. 2011;Corral-Santana et al. 2016).In the interim, XRTs hibernate for decades in a quiescent state, with typical luminosities LX ≲ 10 32 erg s −1 .Optical/near-infrared studies performed during quiescence typically reveal the spectrum of the companion star, which is used as a probe to weigh the accreting star (see Casares & Jonker 2014 for a review of the techniques and inferred BH masses).
The latest version of the BlackCAT catalogue ⋆ E-mail: jorge.casares@iac.es(https://www.astro.puc.cl/BlackCAT/index.php;see Corral-Santana et al. 2016) contains 70 BH XRTs discovered since the start of X-ray missions, with 19 BHs dynamically confirmed (i.e. with dynamical masses in excess of 3 M⊙, the upper bound of a non-rotating neutron star).In addition, by using a scaling relation between the full-widthat-half-maximum (FWHM) of the Hα line (formed in the accretion disc) and the radial velocity semi-amplitude of the companion star K2 (see Casares 2015) indirect evidence for dynamical BHs has been reported in four other XRTs: Swift J1357.2-0933(Mata Sánchez et al. 2015), KY TrA (Zurita et al. 2015), Swift J1753.5-0127(Shaw et al. 2016) and.MAXI J1659-152 (Torres et al. 2021).In all these cases the donor stars are undetected because the optical counterparts are either too faint (r ≳ 23) or strongly contaminated by the light of the accretion disc or a close interloper.Other properties of the Hα line can also be exploited to obtain reliable estimates of the donor to compact star mass ratio q = M2/M1 and the binary inclination (Casares 2016;Casares et al. 2022).In the current paper, we apply this same strategy to the historical XRT GRS 1716-249 (=N Oph 93), which has a challenging quiescent magnitude r=23.2.
GRS 1716-249 (GRS1716 hereafter) was discovered in September 1993 during an X-ray outburst with maximum flux FX (20-100 keV) ∼1.4 Crab (Ballet et al. 1993;Harmon et al. 1993).It was promptly classified as a strong BH candidate on the basis of its X-ray and radio properties, reminiscent of those exhibited by dynamical BHs in the hard state.In particular, GRS1716 displayed a hard power-law spectrum, strong variability at short time-scales, the presence of a lowfrequency type-C quasi-periodic oscillation (QPO), a flat radio spectrum and bright radio flares, correlated with X-rays (Della Valle et al. 1994;van der Hooft et al. 1996;Hjellming et al. 1996;Revnivtsev et al. 1998).The optical counterpart (N Oph 93, V2293 Oph) was identified with a V∼ 16.6 star two days after the outburst peak.No quiescent counterpart brighter than V ∼ 21 could, however, be detected in pre-outburst UK and European Southern Observatory (ESO) Schmidt plates (Della Valle et al. 1994).A tentative ∼14.7 h superhump1 period was proposed from sparse photometry collected in ∼10-60 min V -band nightly blocks over ten consecutive nights during outburst (Masetti et al. 1996, M96 hereafter).The same authors proposed a BH mass ≥ 4.9 M⊙, while Della Valle et al. (1994) argue for a distance of 2.4±0.4 kpc after comparing the peak optical brightness with other BH XRTs.Furthermore, GRS1716 exhibited renewed activity in 1994 and 1995 at 10-20 per cent level of the 1993 main outburst, before settling down in a quiescent state (Churazov et al. 1994;Borozdin et al. 1994;Harmon et al. 1994;Borozdin et al. 1995).
After dwelling for 23 years in quiescence, GRS1716 was found to be active again in December 2016 (Masumitsu et al. 2016;Negoro et al. 2016).Swift and NuSTAR spectral analysis of the new outburst showed this time evidence for a weak thermal disc component.A horizontal displacement across the hardness-intensity diagram towards softer colours was also witnessed, although failing to make a full transition to the canonical soft state (Bassi et al. 2019;Jiang et al. 2020).Type-C QPOs were detected once more in the hard state (Bharali et al. 2019), a classic signature of BH XRTs.Meanwhile, simultaneous radio observations found the system in the radio-quiet branch of the radio/X-ray luminosity plane through the entire outburst, following LR ∝ LX 1.4 (Bassi et al. 2019).On the other hand, optical spectroscopy obtained over the course of the outburst decay provided evidence for conspicuous disc outflows and a possible inflow signature, perhaps caused by a failed wind (Cúneo et al. 2020).Finally, in September 2017 GRS1716 returned to quiescence, where it has remained ever since (Bassi et al. 2019).
Recently, Saikia et al. (2022) reported a detection of the quiescent counterpart of N. Oph 93 with an i ′ -band magnitude of 21.39±0.15.Based on the global optical and X-ray properties observed during the 2016-2017 outburst they also challenge the 2.4 kpc distance of Della Valle et al. (1994) and propose a new value in the range ∼4-8 kpc.Here in this paper we present photometric evidence of a 6.7 h orbital period, based on superhump light curves, recovered from the 1995 mini-outburst, plus an ellipsoidal modulation obtained through deep i-band photometry in 2021.In addition, we report on the presence of an Hα emission line in quiescent spectroscopy and use it to derive dynamical constraints on GRS1716.Our photometric and spectroscopic data strongly advocate for the presence of a BH in GRS1716, and lend support to a binary distance of ∼6.9 kpc, in good agreement with the results of Saikia et al. (2022).

Photometry
In an attempt to measure the orbital period of GRS1716, we first rescued r-band photometry obtained by us with the Danish 1.54-m telescope at the La Silla Observatory in Chile during the 1995 mini-outburst.Observations were taken with the Tektronix TK1024M CCD on the nights of May 7-9 1995.A total of 147×300 s images were collected with the Gunn Rband filter over the three nights, with an average coverage of 5 h per night.Nights were clear, with a median seeing varying between 0.8 and 1.2 arcsec.Differential photometry was performed relative to the local comparison star Pan-STARRS ID 77972599165136381, which has similar brightness to the target.Unfortunately, the reduced images are lost and the raw data can no longer be retrieved since ESO does not keep an archive for the Danish telescope.However, we have managed to recover the data points after digitizing our printed plots of the three nightly light curves.Based on the magnitude of the comparison star quoted by Pan-STARRS (Tonry et al. 2012) we estimate that GRS176 had a mean R-Johnson magnitude of 17.05±0.02at the time.
Photometric observations of GRS1716 were also programmed during the quiescent phase in 2019-2021.First, we obtained exploratory images using the Auxiliary-port CAMera (ACAM) on the 4.2 m William Herschel Telescope (WHT) at the Roque de los Muchachos Observatory in La Palma, on the nights of July 24 and Aug 23 2019.ACAM has a circular field of view with 8.3 arcmin diameter and it is equipped with a low-fringing 2k×4k EEV CCD that gives a platescale of 0.25 arcsec pix −1 .One 600 s integration in Sloan r and another one in Sloan i were obtained on July 24 and three 300 s i-band images on Aug 23under photometric conditions, with seeing ≃0.9 arcsec.
Subsequently, an intensive photometric campaign was performed with the Goodman High Throughput Spectrograph in imaging mode on the Southern Astrophysical Research (SOAR) 4.1m Telescope at Cerro Pachón in Chile.A total of 82×600 s integrations were obtained with the Sloan-i filter on the nights of May 31 and June 1 2021, covering ∼9 h per night.One Sloan-r image was also collected each night to measure the colour of GRS 1716.We employed the red camera of the instrument, a 4k×4k e2v 231-84 CCD which provides a circular field of view of 7.2 arcmin diameter.We used 2×2 binning which results in a platescale of 0.3 arcsec pix −1 .The SOAR observations had to be interrupted for about 100 min every night due to telescope pointing limitations at altitude > 88 • caused by its Alt-Az design.The first frames of each night were affected by problems with telescope tracking and the active optics system.These were rejected as they proved to be unusable for aperture or point spread function (PSF) photometry because of the poor image quality and crowded nature of the field close to the target.This leaves a total of 44 useful images for the photometric analysis with 0.9-1.2arcsec seeing.A log of the photometric observations is presented in Table 1.
The ACAM and SOAR images were reduced in the standard way, with bias subtraction, flat-field correction and stellar alignment using IRAF2 tasks.For the flat-field correction we used a set of sky images obtained during evening twilight.We extracted the stellar counts using PSF photometry with DAOPHOT (Stetson 1987).More details of this process are given in Sect.3.2.2.

Spectroscopy
We also collected spectroscopic observations of GRS1716 on the night of June 22 2020 using the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS; Cepa et al. 2000) at the 10.4m Gran Telescopio Canarias (GTC).These observations were obtained as part of the DDT Director's Discretionary Time program GTC2020-153.Four 1800 s spectra were collected using grism R1000R and a 1.0 arcsec slit, resulting in a wavelength coverage of 5000-10000 Å at 8.3 Å resolution (378 km s −1 at Hα).The slit was oriented at parallactic angle to minimize light losses across the spectrum caused by atmospheric refraction.Seeing conditions were ≤0.9 arcsec through the course of the integrations.The 1-D spectra were extracted using optimal extraction routines within PAMELA (Marsh 1989) after standard debiasing and flat-field corrections.Observations of HgAr+Ne lamps were employed to derive a pixelto-wavelength calibration through a 4th order polynomial fit to 37 lines across the entire wavelength range.Small flexure corrections, obtained from the position of the O i 5577.34 and 6300.30Å sky lines, were applied to individual spectra.The resulting spectra were imported to MOLLY, resampled to a uniform velocity scale of 104 km s −1 and shifted into the heliocentric velocity frame.Each spectrum was also rectified through division by a 3rd order spline fit to the continuum level.

A superhump modulation in 1995
Fig. 1 presents the digitized light curve obtained during the 1995 mini-outburst, where a modulation with a timescale of a few hours is clearly visible.We have computed a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) of these data in the frequency range 0.1 to 10 cycles d −1 at a resolution of 0.001 cycles d −1 and the result is plotted in the lower panel of Fig. 1.The highest peak is found at 3.584 cycles d −1 or 6.7 h and it is marked by a vertical dotted line.The other main peaks correspond to the 1-day alias pattern introduced by our observing window.Note that the 14.7 h period claimed by M96 (indicated in the figure with an arrow) has significantly lower power as it is only the 4th highest peak in the periodogram.Inspection of the light curve, folded on the four main periods (Fig. 2 ) strongly favours 0.279 d (=6.7 h) as the true period since it is the one that shows minimum dispersion.A Gaussian fit to the highest peak in the Lomb-Scargle periodogram yields 0.279±0.009d, where we have adopted the sigma of the Gaussian as a conservative estimate of the statistical uncertainty.
Regarding the discrepancy with M96's period, we note that the latter is based on eight V -band night visits covering ∼20-60 min each (plus two additional 10-15 min visits), which provides a rather poor phase coverage.The 14.7 h period is tentatively selected by M96 after substantial data processing through iterative cleaning of the frequencies introduced by the window function.Conversely, the 6.7 h period is almost completely sampled on two of our three consecutive visits and stands out as the peak with the highest power in the dirty periodogram, without any need for data processing (Fig. 1).Furthermore, the 14.7 h period is clearly disproved by the large scatter shown by our data in the folded light curve.All in all we conclude that M96's period is an alias of the 6.7 h modulation and should not be used in future studies of GRS1716.
The light curve folded on the 0.279 d period exhibits a sawtooth morphology, with a small secondary peak around minimum light that moves in phase across nights.Triangular modulations, with superimposed (higher frequency) waves are typically observed in outburst light curves of BH XRTs and are interpreted as superhumps (O'Donoghue & Charles 1996;Zurita et al. 2008;Thomas et al. 2022).As a matter of fact, the extreme mass ratios, characteristic of BH XRTs, tend to prevent donor stars from being exposed to X-ray irradiation and, consequently, to exhibit orbitally modulated light curves (see e.g.Torres et al. 2021).In any case, since superhump periods are typically ≈2 per cent longer than orbital, the latter (if present) would not be resolved from the superhump frequency in our data.

WHT+ACAM 2019 Epoch
Fig. 3 (left panel) displays the ACAM Sloan-r band image of the GRS1716 field, obtained on the night of July 24 2019.The field was calibrated relative to a reference star (labelled R in the figure), selected from the Pan-STARSS Data Release 2 catalogue (Chambers et al. 2016).Its magnitude and colour are listed in Table 2. Due to the faintness of the target and the presence of nearby bright stars we used point-spread function (PSF) photometry to measure the brightness of GRS1716 and find r = 22.60 ± 0.06.A consecutive i-band image, collected only 10 min after, gives i = 21.59 ± 0.03.On the other hand, the mean magnitude of three i-band images obtained on Aug 23 2019 gives i = 22.08 ± 0.07, i.e. a 0.5 mag difference which may be ascribed to a combination of flickering, orbital phasing and possible episodic variability between the two visits.

SOAR+Goodman 2021 Epoch
The right panel of Fig. 3 shows a close-up view of the GRS1716 field in the i-band, observed with SOAR on 2021 May 31.Standard PSF photometry was also performed during the SOAR campaign.A Moffat distribution model was constructed for seven selected stars after removing the neighbours.This was fitted to the target and a number of close stars in the field, with their positions fixed to avoid the PSF process shifting the position of faint stars towards brighter nearby objects.After several tests, we set the fitting and aperture radii to 2 pixels.Differential photometry was performed on the target and two comparison stars (labelled as C1 and C2 in Fig. 3), relative to a local standard (S) whose Pan-STARRS magnitude is i = 20.68 ± 0.04.Note that stars brighter than S (e.g.reference star R) were registered outside the detector linearity regime.Comparison C2 has similar brightness to the target while C1 was chosen to monitor the  possible impact that the bright star, located 2.1" SW of the target, could have in the PSF photometry.Table 2 summarizes the mean magnitudes of all the reference stars.
Fig. 4 shows the light curve of GRS1716 and stars C1 and C2.Unlike the comparison stars, GRS1716 is clearly variable.As is usually the case in quiescent XRTs, we expect this variability to be dominated by the changing visibility of the tidally distorted companion star (the so-called ellipsoidal modulation), with some contamination due to the accretion disc light (Martin et al. 1995;Shahbaz et al. 1996).The ellipsoidal modulation has two maxima and two minima per orbital cycle and its amplitude increases with inclination.To determine the orbital period of GRS1716 we have computed a Lomb-Scargle period analysis of the SOAR light curve and the result is displayed in the bottom panel of Fig. 4. Note that, because of the double modulation per orbital cycle of the ellipsoidal variability we opt for plotting half the frequency in the x-axis.
The highest peak (0.324 d) is likely an artefact of data sampling because the folded light curve only covers ∼60 per cent of the orbit (see top panel in Fig. 5).In addition, it does not match any of the strong peaks in the 1995 outburst periodogram (Fig. 1).As a matter of fact, the second highest peak in the SOAR periodogram coincides with the superhump frequency (ν/2=3.597cycles day −1 or 0.278 d), a strong indication that this is the true orbital period (P orb ).The remaining peaks at shorter frequencies also show signifi- cant gaps in their phase folded light curves and we interpret them as sidelobes of the orbital frequency.Further independent support for the 0.278 d orbital period is provided by the detection of a ≈6.5 h modulation in both the radial velocities and the EW of the Hα line observed during the 2017 outburst (see Appendix A).We hereafter adopt 0.278 ± 0.008 d as the orbital period in GRS1716, where the error corresponds again to the sigma of a Gaussian fit to the peak in the periodogram.
The bottom panel in Fig. 5 displays the GRS1716 quiescent light curve folded on P orb =0.278 d.We have adopted a tentative zero phase T0(HJD)=2459365.963which places the deepest minimum at orbital phase 0.5, referred to as the superior conjunction of the companion star.The reason for a deeper minimum lies in the larger gravity darkening of the inner Lagrangian point.The large amplitude of the folded light curve and the difference between the two minima suggests that GRS1716 is seen at high inclination.In order to constrain the latter we have performed an ellipsoidal model fit using the code XRBINARY (Bayless et al. 2010).This is a staronly fit, where we have assumed a binary mass ratio q = 0.1, K2 = 521 km s −1 (Sect.3.3) and T eff = 4300 K (i.e.adequate for a K6V companion, see Sect.4.1.1).The flux spectrum of the companion star is computed from the stellar atmosphere models of Kurucz (1996), using a non-linear limb-darkening law (Claret 2000b) valid for log g = 0.0 − 0.5.The gravity darkening depends on the star's effective temperature and is based on Claret (2000a).The best fit is obtained for i = 72 +4 −5 deg, although this number should be treated with caution.To start with, the ellipsoidal modulation appears strongly contaminated by flickering activity, as is commonly observed in the light curves of quiescent XRTs (e.g.Zurita et al. 2003).In addition, its amplitude might be diluted by the (steady) contribution of the disc light, which has so far been neglected.
In an effort to estimate the impact of flickering in our modeling we have followed Pavlenko et al. (1996) and fitted the lower envelope of the light curve.This is expected to trace the true ellipsoidal variability, with minimum flickering contamination.We computed the lower envelope curve through an iterative sigma clipping process that rejects data points > 2σ above a double sine wave model, with periods fixed to P orb and 0.5×P orb .The latter is a good approximation to the real ellipsoidal modulation, but much less expensive in terms of computing time.The process converges after four iterations.A proper ellipsoidal fit to the resulting lower envelope curve now yields i = 68 ± 6 • (see bottom panel in Fig. 5).This exercise suggests that flickering variability tends to bias inclinations high while also underestimates the quoted uncertainties.

Quiescent magnitudes and colour
The mean i-band magnitude of the two SOAR nights is i = 22.14 ± 0.14, where the error reflects the standard deviation of the light curve.Note that this is significantly fainter than the quiescent magnitude reported in Saikia et al. (2022).In order to infer the quiescent colour of GRS1716 we examined the two r-band images gathered in this run and find r = 23.20 ± 0.10 (2021 May 31) and r = 23.70 ± 0.20 (2021 June 1).Given the large variability exhibited by GRS1716 we need to compare these magnitudes with those from the nearest iband images and find (r − i) = 1.10 ± 0.13 and (r − i) = 1.30 ± 0.22 respectively.These values are consistent within errors with the measurement (r − i) = 1.01 ± 0.07, obtained on 2019 July 24 with ACAM.This indicates that, despite intrinsic variability, the colour appears to be very stable over the years.The weighted mean of the three values, i.e. (r−i) = 1.05 ± 0.06, will be adopted hereafter as representative of the quiescent colour of GRS1716.Similarly, we combine this colour with the average i-band magnitude of the two SOAR nights (i = 22.14 ± 0.14) to derive a mean quiescent r-band magnitude of 23.19±0.15.

Quiescent spectroscopy
The four GTC spectra collected in 2020 have a mean signalto-noise ratio of ≈2 in the continuum.The only notable feature in the average spectrum is the presence of a broad double-peaked Hα emission line, characteristic of an accretion disc (Fig. 6).The line is strong, with a mean equivalent width EW=107±5 Å, and its centroid velocity is redshifted by 209 ± 54 km s −1 .Since no absorption lines from the companion are detected we need to exploit the properties of the Hα profile to set constraints on the system parameters.
First, we use the FWHM-K2 correlation to determine the radial velocity semi-amplitude of the donor star K2 (Casares 2015).The width of the line is measured by fitting a single Gaussian to the four individual spectra in a window of ±10000 km s −1 centered on Hα.This yields an average of FWHM=2149±260 km s −1 , where the quoted errors repre- sent the standard deviation in the distribution of values.Since the individual spectra are rather noisy we have also fitted the averaged spectrum and find FWHM=2235±137 km s −1 .Observations of quiescent XRTs show evidence of accretion disc activity, leading to ∼10 per cent variability in FWHM (Casares 2015).Given the limited phase coverage (∼25 per cent of the 0.278 d period) we decided to adopt the results of the fit to the averaged spectrum, but with a conservative 10 per cent error to account for possible orbital and episodic variability i.e.FWHM=2235±224 km s −1 .The empirical scaling K2 = 0.233 × FWHM thus yields K2 = 521 ± 52 km s −1 .
Second, we use the correlation between q and the ratio of the double peak separation (DP ) to the line width (FWHM) to estimate the former (Casares 2016).A symmetric 2-Gaussian model fit to the average profile gives DP = 1263 ± 49 km s −1 , which leads to q = 0.074 +0.349  −0.059 (68 per cent confidence through a Monte Carlo simulation, see Casares 2016) via the equation log q = −6.88− 23.2 log(DP/FWHM).Unfortunately, the large errors caused by the limited quality of the Hα profile leaves this parameter largely unconstrained.
Finally, we measure the depth of the Hα central trough to constrain the inclination as in Casares et al. (2022).The same 2-Gaussian model yields W = 973 ± 72 km s −1 for the full-width-half-maximum of the two fitted Gaussians which, together with DP , leads to a normalized double peak trough depth T = 1 − 2 1−(DP/W ) 2 = 0.40 ± 0.12 and thus i = 61 ± 12 • via equation i(deg) = 23.7 + 93.5 T .Here, T has been corrected by ∆T = 0.097 × ∆V 2 res /(DP 2 − W 2 ) = +0.02 in order to account for smearing in the Hα profile caused by our limited spectral resolution ∆Vres = 378 km s −1 .The quoted uncertainty in the inclination has been computed through a Monte Carlo simulation as in Casares et al. (2022).

DISCUSSION
By combining our photometric determination of the orbital period with the K2 velocity of the companion star we constrain the mass function of the compact object to This provides strong support for the presence of a BH as it exceeds the ≈3 M⊙ limit for a stable neutron star.The BH nature is consistent with the global X-ray properties, in particular the detection of Type-C QPOs in both the 1993 and 2017 outbursts (van der Hooft et al. 1996;Bharali et al. 2019).Note that our mass function relies on a tight FWHM-K2 correlation that stems from a fundamental scaling between the dynamics of the accretion disc and the velocity of the companion star in quiescent XRTs (Casares 2015).Arguably, this is the only avenue to derive dynamical information on such a faint (r ≈23) target, as traditional methods (i.e.Doppler monitoring of weak metallic lines from the companion) are exceedingly time consuming or simply impossible with the available instrumentation.For a better knowledge of the BH mass we now need to look at other observational constraints that can be enforced on the binary parameters.

Spectral type of donor star
Under the assumption that the companion star fills its Roche lobe, the orbital period implies a mean stellar density ⟨ρ⟩ ≈ 2.5 gr cm −3 through equation ⟨ρ⟩ ≈ 110 × P −2 orb , where P orb is in units of hours (Frank et al. 2002).Such density is consistent with a ≈K5 main sequence star (Drilling & Landolt 2002).The spectral type of the donor star can, in principle, be constrained through the observed quiescent (r − i) colour.Our photometry gives (r − i) = 1.05 ± 0.06 which, after correcting for interstellar (IS) extinction in the Sloan bandpasses (Ar = 2.285 × E(B − V ) and Ai = 1.698 × E(B − V ); see Schlafly & Finkbeiner 2011), together with our reddening determination E(B − V ) = 1.0 ± 0.1 (see Appendix B), results in (r − i)0 = 0.46 ± 0.29.This agrees well with a K6 main sequence star (Covey et al. 2007) although, strictly speaking, it only sets an upper limit to the spectral type since we have so far neglected any (blue) light contamination from the accretion disc.For example, Cantrell et al. (2010) estimate the fractional contribution of the accretion disc to the total light in A0620-00 during the passive-state (i.e. that with the lowest flux level and variability), both in the V and I bands as fV = 0.35 and fI = 0.25, respectively.A linear interpolation to the effective wavelength of the Sloan filters yields fr = 0.32 and fi = 0.27, which leads to excess magnitudes ∆r = 0.42 and ∆i = 0.34.If our GRS1716 magnitudes were affected by a similar amount of disc contamination as in the passive-state of A0620-00, the dereddened (r−i) colour would rise by ∼0.08 mag, thus favouring a ∼K7 spectral type.It should be borne in mind that quiescent BH XRTs are sometimes caught in the so-called active-state, with higher levels of disc activity, although we note that this does not seem to have a large impact on the observed colour (Cantrell et al. 2008;Dinçer et al. 2018).In conclusion, after correcting for extinction and disc veiling, we find that the quiescent (r − i) colour of GRS1716 favours a ≈K6-7 spectral type, even though we caution that large uncertainty is associated with the knowledge of IS reddening.Incidentally, this classification is in excellent agreement with the spectral type predicted by the empirical spectral type-period relationship of Smith & Dhillon (1998) for P orb =6.7 h.

Mass ratio
The evidence for superhumps in 1995 implies that the accretion disc expanded beyond the tidal 3:1 resonance radius and, therefore, q ≲0.25 (Paczyński 1977).An independent constraint is placed by the spectral classification of the companion star.Evolutionary models of short period (≲12 h) Xray binaries provide an upper limit to the mass of the donor star based on zero-age main sequence stars with the same T eff (Kolb et al. 2001).Our K6 V classification thus implies M2 ≲ 0.69 M⊙ and, since the compact object is a BH with M1 ≳ 3 M⊙, then q ≲ 0.23.A tentative mass ratio can in fact be derived from the radial velocities of the Hα line during outburst.We show in Appendix A that the red peak of the line on 2017 June 3 traces a sine wave with a semi-amplitude of 79 km s −1 .Under the assumption that this reflects the motion of the compact star we derive q ≈ 79/521 = 0.15, although we warn about possible biases in the radial velocity curve caused by contamination from a putative superhump (e.g.Zurita et al. 2002;Torres et al. 2002;Haswell et al. 2004).Finally, our 2-Gaussian fit to the quiescent GTC spectrum gives a very loose constraint q = 0.015 − 0.423 (68 per cent confidence).It should be noted that large q values, such as 0.423, can in fact be ruled out because, when combined with M1 ≳ 3 M⊙, imply M2 ≳ 1.3 M⊙ i.e. a F5 V star that would not fit in the Roche lobe of a 6.7 h orbit.Furthermore, the ensemble of q values for BH XRTs with low-mass companions are narrowly distributed around q = 0.06, and none exceeds q = 0.18 (Casares 2016;Heida et al. 2017;Torres et al. 2020).To summarize on the mass ratio restrictions, given the information that is available from the photometric evidence of superhumps, the orbital period, the Hα velocities and double peak fit we decide to adopt the wide range q = 0.02 − 0.25.

Inclination
The 2017 outburst of GRS1716 was intensively scrutinized by Swift and NuSTAR, with no X-ray eclipses nor dips being detected (Bharali et al. 2019).This sets a rough upper limit to the binary inclination i ≲ 70 • (Frank et al. 1987).On the other hand, the observation of double-peaked Hα profiles in quiescence and outburst suggests i ≳ 50 • (see Fig 6 and Appendix A).Further, by fitting the depth of the double peak trough in the quiescent GTC profile we measure i = 61 ± 12 • (Sect.3.3).This agrees well with the inclination obtained for XTE J1859+226 (i = 66.3 ± 4.3 • ), a twin of GRS1716 with identical P orb and Hα line properties (Corral-Santana et al. 2011;Yanes-Rizo et al. 2022).In particular, both lines are equally strong (EW=102±12 Å and 107 ± 11 Å, respectively), a fact that suggests small foreshortening of the disc continuum light and, therefore, a moderately high inclination in the two systems (Warner 1986).A relatively high inclination i = 72 +4 −5 deg is favoured by the ellipsoidal fit to the quiescent i-band light curve, although systematics caused by flickering might be important.Circumstantial evidence for a relatively high inclination is also provided by the large disc wind velocities observed by Cúneo et al. (2020) during the 2016-2017 outburst (see Table 2 in Panizo-Espinar et al. 2022).Given all the above we decide to take the conservative value i = 61 ± 15 • , where we have adopted a larger uncertainty than provided by the fit to the Hα trough in order to encompass the inclination range allowed by the ellipsoidal modeling (Sect.3.2.2).

Black hole mass
In order to estimate the BH mass in GRS1716 we have run a Monte Carlo simulation with 10 5 trials on eq. 1 and the constraints above.In summary, we take Gaussian distributions for K2 and P orb with mean and standard deviations set by our observed values i.e.K2 = 521 ± 52 km s −1 and P orb = 0.278±0.008d.We also adopt a normal distribution of inclination angles with i = 61 ± 15 • , but truncated at i = 90 • .Regarding the mass ratio we assume a uniform distribution between q = 0.02 − 0.25.This yields M1 = 8.3 +7.4 −3.1 M⊙ and M2 = 1.1 +1.4  −0.7 M⊙ for the masses of the two binary components.We note, however, that most M2 solutions are ruled out by our spectral type determination and, therefore, decided to enforce further ad hoc restrictions.First, from the lower limit of our quiescent de-reddened colour (r − i)0 = 0.46 ± 0.29 we set a conservative upper limit to the donor's spectral type of G4V and thus M2 < 0.98 M⊙.Subsequently, given the absence of X-ray eclipses we impose a more restrictive upper limit on inclination based on the geometrical limit cos i < R2/a, where the ratio between the donor star radius R2 and binary separation a is tied to q by Eggleton's expression (Eggleton 1983).With these additional constraints we find M1 = 6.4 +3.2 −2.0 M⊙ and M2 = 0.5 ± 0.3 M⊙ through the Monte Carlo analysis.
The large fractional uncertainty in K2 and, especially, the loose constraint on inclination lead to a BH mass distribution that is very skewed towards high values.Tighter constraints on these two parameters are thus fundamental to obtain a more accurate BH mass.An independent determination of the BH mass can also be obtained by simply combining the intrinsic width of the double peak profile (W ) with P orb through equation 8 in Casares et al. ( 2022) where β = 0.84 accounts for the fact that the outer disc material is ≈84 per cent sub-Keplerian.This yields M1 = 7.0 ± 1.5 M⊙, where, following Casares et al. (2022), we adopt a 20 per cent error in M1, based on the distribution of differences with respect to the dynamical masses of a sample of calibrators.

Distance, Galactic elevation and X-ray luminosity
Our quiescent r-band magnitude and orbital period allow us to constrain the luminosity of GRS1716 and its distance.In order to do so we apply the empirical Mr − P orb correlation obtained for quiescent BH XRTs in Casares (2018) (eq.8): Mr = 4.64(0.10)− 3.69(0.16)log P orb (d), For P orb = 0.278 d we get Mr = 6.7 ± 0.2 which, in turn, is fully consistent with the absolute r-band magnitude of a K5-7 main sequence star.The distance to GRS1716 is obtained by comparing Mr with the observed magnitude r = 23.19 ± 0.15 through the usual distance modulus equation We here adopt an r-band extinction Ar = 2.285 E(B − V ) (Schlafly & Finkbeiner 2011), which assumes a traditional reddening law RV = AV/E(B − V ) = 3.1, appropriate for the diffuse interstellar medium (Cardelli, Clayton & Mathis 1989;Fitzpatrick 1999).From the strength of IS lines in the 2017 outburst spectra we find E(B − V ) = 1.0 ± 0.1 (see Appendix B) which results in Ar = 2.29 ± 0.23 and a distance d = 6.9 ± 1.1 kpc.We note that this is a factor ≈ 3 larger than the 2.4±0.4 kpc reported in Della Valle et al. (1994), but agrees well with the distance favoured by Saikia et al. (2022).
For an independent test on our distance we have looked at 3D reddening maps constructed from Pan-STARRS/2MASS photometry and GAIA parallaxes (Green et al. 2019).These show a steady rise in extinction up to 2.5-3.0 kpc, where it jumps to Ar = 2.617 E(g − r) = 2.4.Beyond this point, the extinction level remains constant all the way up to ≈6.5 kpc.We therefore conclude that our estimated distance and reddening values are broadly consistent with the available 3D reddening maps although unfortunately the latter are not very informative beyond ≳3 kpc.
Our distance determination has wide-ranging implications.To start with, it places the binary at z ∼0.8 kpc above the Galactic Plane, in better agreement with an observed P orb −z anticorrelation which lends support to the existence of natal kicks (Gandhi et al. 2020).In relation to this, Atri et al. (2019) estimated a potential kick velocity (i.e. the peculiar velocity when the system crosses the Galactic Plane) for GRS1716 of 67±21 km s −1 .This was based on its VLBI proper motion, vertical elevation and the Galactic rotation velocity for an assumed distance of 2.4 kpc.Our revised distance and z-elevation, together with a tentative systemic velocity of 209 ± 54 km s −1 (derived from the centroid of the quiescent Hα line), implies a larger Galactic Plane crossing velocity in the range ∼153-355 km s −1 (90 per cent confidence; P. Atri, private communication), lending further support for a strong natal kick in GRS1716.Next, the new distance raises the peak X-ray and radio luminosities by a factor ∼8.3.Although these still keep GRS1716 on the radio quiet branch of the LR-LX diagram (Bassi et al. 2019;Saikia et al. 2022), the new peak X-ray luminosity during the 2016-2017 outburst becomes LX(0.5-10 keV)=4.4×10 37ergs s −1 .This corresponds to ≈ 6% L Edd (for a 6 M⊙ BH), in good agreement with the linear P orb − LX correlation proposed by Wu et al. (2010).Such luminosity is also more consistent with the presence of a thermal disc component, detected during the hard-intermediate state, and the evidence that GRS1716 attempted a failed transition to the soft-state (Armas Padilla & Muñoz-Darias 2017; Bassi et al. 2019;Jiang et al. 2020).In fact, hard-to-soft state transitions in BH XRTs are typically observed at ≈0.1 L Edd but rarely at ∼ 10 −3 L Edd (Maccarone 2003;Yu & Yan 2009;Dunn et al. 2010), as would be implied if the distance were 2.4 kpc.A larger distance also shifts GRS1716 from the neutron star to the BH region in the optical/X-ray correlation plots (see Saikia et al. 2022).And finally, the new distance will surely affect the accretion disc parameters derived through X-ray spectral modeling, such as the mass accretion rate, disc inclination and inner disc radius.For example, it increases the mass accretion rate ṁ (=L 0.1−100kev /0.2 L Edd ) by one order of magnitude, thus helping to reconcile the low disc density measured by Jiang et al. (2020) with that of GX 339-4 (at similar ṁ).Our restrictions will certainly have an impact too on the BH spin inferred from spectral fitting (Tao et al. 2019).Clearly new analysis of the 2017 X-ray spectra should be performed in the light of the new BH mass and distance derived in this paper.

CONCLUSIONS
We have presented evidence for a 0.278±0.008d orbital period in the BH XRT GRS 1716-249, based on 1995 superhump variability and ellipsoidal light curves obtained in quiescence.The extreme faintness of the optical counterpart (r = 23.19 ± 0.15) prevents a direct detection of spectral features from the companion star in quiescence.However, GTC spectroscopy reveals a strong double-peaked Hα emission line that we exploit to derive fundamental binary parameters through a set of scaling relations (Casares 2015(Casares , 2016;;Casares et al. 2022).Using these, combined with light curve information, leads to f (M ) = 4.1 ± 1.2 M⊙, q = 0.02 − 0.25 and i = 61 ± 15 • .By including further restrictions on the companion's spectral type and the lack of X-ray eclipses, we constrain the BH mass to M1 = 6.4 +3.2 −2.0 M⊙.Finally, we find that GRS 1716-249 is located 6.9±1.1 kpc away and 0.8 kpc above the Galactic Plane, in support for a strong natal kick.Our new distance implies that the 2017 failed transition to the soft state occurred at ≈ 6% L Edd , in agreement with observations of other BH XRTs.2020) for technical details on the acquisition and reduction of the data.All the spectra have been normalized to the continuum level.The May 19 spectrum displays a classic double-peaked profile , characteristic of an accretion disc (Fig. A1).The line is moderately narrow, with FWHM=1106±9 km s −1 , as measured through a Gaussian fit.We also derive a double peak separation of 650±3 km s −1 after fitting a double Gaussian model to the profile.
While Cúneo et al. (2020) focus on the detection of outflow signatures in the Balmer lines we here search for evidence of the orbital period by looking at variability in the radial velocities and EW of the Hα emission, the strongest line in the entire spectrum.The ten 916 sec spectra from June 3 span over a continuous interval of ∼6 h and therefore nearly cover a full orbital cycle.Fig. A2 presents the trailed spectra of the Hα line on that night, with two orbital cycles being plotted.We note that the blue part of the line is strongly affected by transient P-Cygni absorptions (see the skewed profiles in Fig. 4 of Cúneo et al. 2020) although the red peak remains largely unaffected.Visual inspection shows that the latter traces a clear sinusoidal wave.
In order to extract the red peak velocities we performed a two-Gaussian model fit to the individual profiles.The model consists of a broad base and a narrow redshifted component with all the parameters left free.The EW of the entire profile was also measured by integrating the line in a band of ±1500 km s −1 around the rest velocity.Fig. A3 displays the red peak velocities and EWs of the Hα line, folded on the 0.278 d orbital period.Although our limited time coverage prevents us from measuring an accurate period, the observed variability is clearly consistent with 0.278 d and not with the 0.613 d period reported by M96, unless the Hα velocity semiamplitude were implausibly large at ∼235 km s −1 .Sine wave fits to the red peak velocities and EWs, with the period fixed to 0.278 d, gives semi-amplitudes of 79±2 kms and 1.49±0.02Å, respectively.

APPENDIX B: INTERSTELLAR EXTINCTION
To constrain the reddening towards GRS1716 we have measured the strength of several IS absorption lines in the average of all the X-shooter spectra.For the case of the IS line at λ6613 (Fig. A1) we find EW=0.20±0.01Å which, according to the linear calibration of Herbig (1975), corresponds to E(B − V ) = 0.8 ± 0.2.The quoted error in EW is a conservative estimate based on a range of possible continuum positions, while the uncertainty in E(B − V ) is determined by the range of fitting coefficients listed in Table 4 of Herbig (1975).On the other hand, the Na D1 line at λ5896 yields EW=0.64±0.1 Å which translates into E(B −V ) > 0.75 (Munari & Zwitter 1997).We take this value as a conservative lower limit because the line is most likely saturated.Evidence for this is provided by the ratio EW(D2)/EW(D1)=1.1.This should be ∼2 at low optical depths while it converges asymptotically towards 1.1 at large reddening values.
The K I line at λ7699 offers a more robust diagnostic for the case of high reddening values (Munari & Zwitter 1997).We measure EW(K I)=0.26±0.01Å in our spectrum which, according to the empirical fit of Munari & Zwitter (1997) results in E(B − V ) = 1.05 ± 0.05.This value agrees well with the estimate of Della Valle et al. (1994), who derive E(B−V ) = 0.9±0.2 from a combination of different methods, including the strength of several IS lines, the outburst (B − V ) colour, the Balmer decrement and the column density of atomic hydrogen and molecular CO.It also agrees with the hydrogen column density NH = 0.70 × 10 22 cm −2 measured during the 2017 outburst (Bassi et al. 2019(Bassi et al. , 2020)), which implies E(B − V ) = 1.03 (Güver & Özel 2009).Considering all this information, we decide to adopt E(B − V ) = 1.0 ± 0.1 as our best estimate for the interstellar reddening towards GRS1716.

Figure 1 .
Figure 1.Top panel: Light curve of GRS1716 obtained on the nights of 7-9 May 1995.A sine wave with a 0.279 d (=6.7 h) period is overplotted.Bottom panel: Lomb-Scargle periodogram of the 1995 photometry.The red dashed line indicates the highest peak, corresponding to 3.584 cycles d −1 (=6.7 h).The vertical arrow labels the 14.7 h period reported by M96.

Figure 2 .
Figure 2. Light curve of GRS1716 folded on the four periods with the highest power in the Lomb-Scargle periodogram.The 0.279 d period (panel a) is clearly favoured by the reduced data dispersion.Panel d displays the data folded on M96's period.

Figure 3 .
Figure 3. Left panel: 600 s r-band image of GRS1716 in quiescence, obtained with ACAM on the night of 2019 July 24.The target is indicated by a cross mark, while reference star R with an arrow.North is up and East to the left.Right panel: Zoom in on the average of the 10×600 s i-band images of GRS1716 with seeing <1 arcsec, obtained with SOAR on the night of 2021 May 31.The target, two comparison stars (C1 and C2) and a local standard (S) are indicated.

Figure 4 .
Figure 4. Top panel: SOAR light curve of GRS1716 in quiescence and two comparison stars, obtained on the nights of 2021 May 31 and June 1.For the sake of display we have shifted star C2 by -0.6 mag.Bottom panel: Lomb-Scargle periodogram of the GRS1716 light curve.We plot half the frequency in the x-axis because of the double-wave nature of the expected ellipsoidal modulation.The red dashed line indicates the superhump frequency (=6.7 h), which coincides with the second highest peak in the periodogram.

Figure 5 .
Figure 5. Top panel: SOAR light curve of GRS1716 folded on 0.324 d, the period associated with the highest peak in the Fig. 4 periodogram.Bottom panel: same as before, but folded on 0.278 d , our favoured orbital period.The best ellipsoidal model fit with i = 72 • is overplotted in solid line, while an ellipsoidal fit to the lower envelope of the light curve is shown in dashed line style.The latter yields a best inclination i = 68 • (see Sect. 3.2.2 for details).

Figure 6 .
Figure 6.Hα profile of GRS1716 in quiescence, obtained from the average of four 1800 s GTC spectra.A single Gaussian fit is plotted in red dashed-line style and a symmetric 2-Gaussian model in blue dashed line.

Figure A1 .
Figure A1.Hα profile on the night of 2017 May 19.The sharp absorption to the right of Hα corresponds to the λ6613 IS band.

Figure A2 .
Figure A2.Trailed spectra of the Hα line on the night of 2017 June 3.The 10 individual spectra are replicated over two orbital cycles.The line fux has been normalized by its EW to remove large contrast variations, with a factor ∼2 amplitude.These are dominated by changes in the intensity of the red peak.

Figure A3 .
Figure A3.Top: radial velocities of the red peak in the Hα profile.The data are folded onto two cycles, with phase zero arbitrarily assigned to the time of the first spectrum.A sine wave fit with P=0.278 d is overplotted.Bottom: EWs of the Hα line folded onto two cycles, with the best sine wave fit superimposed.Errorbars are always smaller than the symbol size.

Table 1 .
Journal of photometric observations

Table 2 .
Magnitudes of GRS1716 and PS1 reference stars †We quote the mean magnitudes and rms from the nights of 2021 May 31 and June 1.