A Multiwavelength Study of GRS 1716-249 in Outburst: Constraints on Its System Parameters

We present a detailed study of the evolution of the Galactic black hole transient GRS 1716−249 during its 2016–2017 outburst at optical (Las Cumbres Observatory), mid-infrared (Very Large Telescope), near-infrared (Rapid Eye Mount telescope), and ultraviolet (the Neil Gehrels Swift Observatory Ultraviolet/Optical Telescope) wavelengths, along with archival radio and X-ray data. We show that the optical/near-infrared and UV emission of the source mainly originates from a multi-temperature accretion disk, while the mid-infrared and radio emission are dominated by synchrotron emission from a compact jet. The optical/UV flux density is correlated with the X-ray emission when the source is in the hard state, consistent with an X-ray irradiated accretion disk with an additional contribution from the viscous disk during the outburst fade. We find evidence for a weak, but highly variable jet component at mid-infrared wavelengths. We also report the long-term optical light curve of the source and find that the quiescent i′ -band magnitude is 21.39 ± 0.15 mag. Furthermore, we discuss how previous estimates of the system parameters of the source are based on various incorrect assumptions, and so are likely to be inaccurate. By comparing our GRS 1716−249 data set to those of other outbursting black hole X-ray binaries, we find that while GRS 1716−249 shows similar X-ray behavior, it is noticeably optically fainter, if the literature distance of 2.4 kpc is adopted. Using several lines of reasoning, we argue that the source distance is further than previously assumed in the literature, likely within 4–17 kpc, with a most likely range of ∼4–8 kpc.


INTRODUCTION
Black hole X-ray binaries (BHXBs) are interacting binary systems composed of a black hole (BH) accreting matter from a secondary companion star. The accreted matter forms a differentially rotating disk around the BH known as an accretion disk (Shakura & Sunyaev 1973). A large fraction of the accretion energy is often channeled into relativistic, collimated outflows known as jets (e.g. Blandford & Konigl 1979;Fender, Belloni & Gallo 2004). Many BHXBs are transient in nature, alternating between periods of quiescence (typically lasting years to decades, with the X-ray luminosities in the range of 10 30−33 erg s −1 ) and outburst (typically lasting weeks to months, with X-ray luminosities reaching 10 36−39 erg s −1 , e.g. Corral-Santana et al. 2016;Tetarenko et al. 2016).
During an outburst, many BHXBs undergo hysteresis in the spectral state transitions following a q−shaped evolutionary pattern in the hardness-intensity diagram (HID; Miyamoto et al. 1995;Homan et al. 2001;Belloni 2010). The rise of the outburst is generally dominated by a hard, power law-like spectral component (with photon index Γ <2) with an high-energy cut-off at 50-100 keV. This is known as the hard state (HS), which is usually associated with thermal Comptonization due to Compton up-scattering of soft disk photons by a corona of hot electrons (e.g. * NASA Einstein Fellow Thorne & Price 1975;Sunyaev & Titarchuk 1980;Done et al. 2007). During the HS, collimated compact jets are launched, emitting self-absorbed synchrotron emission that dominate radio though infrared (IR) wavelengths (e.g. Corbel et al. 2000;Fender, Belloni & Gallo 2004), in analogy with those observed in active galactic nuclei (Blandford & Konigl 1979;Hjellming & Johnston 1988). Many BHXBs in the HS follow a non-linear radio/X-ray luminosity correlation, where L R ∝ L β X with β ∼0.5-0.7 (e.g. Corbel et al. 2003Corbel et al. , 2013Gallo, Degenaar & van den Eijnden 2018), which extends to active galactic nuclei through the fundamental plane of BH activity (Merloni et al. 2003;Falcke et al. 2004;Saikia et al. 2015Saikia et al. , 2018, suggesting scale invariance of compact jets. During the peak and decay of an outburst, when the system is said to be in the soft state (SS), the spectra are dominated by a soft, blackbody-like spectral component due to an optically thick, geometrically thin accretion disk (Shakura & Sunyaev 1973). The jets are suppressed in this state (e.g. Tananbaum et al. 1972;Fender et al. 1999;Coriat et al. 2011;Russell et al. 2011b;Koljonen et al. 2018;Russell et al. 2019b;Carotenuto et al. 2021). During the transition between these two states, the system enters the intermediate state (IS), dominated by a thermal disk component with a color temperature of 0.1-1 keV, which is further classified based on the Xray timing properties into hard-intermediate and softintermediate states (e.g. Belloni 2010). Depending on the source state, fast variability can be observed, including quasi-periodic oscilla-tions (QPOs), that have been classified into three types: A, B and C (e.g. Ingram & Motta 2019). A number of BHXBs remain in the HS for the entire duration of the outburst (or only transition to the hardintermediate state). These are referred to as 'hard-only state outbursts' (Tetarenko et al. 2016), 'low/hard state outbursts' (Belloni et al. 2002), 'failed outbursts' (e.g. Capitanio et al. 2009;Curran & Chaty 2013), 'failed state transition outbursts' (Bassi et al. 2019) or 'failedtransition outbursts' (Alabarta et al. 2021).

GRS 1716−249
In 1994 September, GRS 1716−249 had a series of several X-ray re-flares or mini-outbursts, as observed by both SIGMA and BATSE at the level of ∼10% of its peak value in 1993 (Revnivtsev et al. 1998). During this period, the X-ray light curve was dominated by at least four sawtooth-like rebrightening events with slow rise (∼30-70 days) and dramatic decay (∼10 days), accompanied by simultaneous radio flares following the onset of decays (Hjellming et al. 1996). This re-brightening event lasted ∼400 days and had at least 4 separate peaks in hard X-rays.
The source had another outburst after almost 21 years in quiescence, and was detected by the Monitor of Allsky X-ray Image (MAXI) on 2016 December 18 (MJD 57740, Negoro et al. 2016), with a photon index of Γ= 1.62±0.06 on 2016 December 21 (MJD 57743, Mashumitsu et al. 2016). It was found to be in the hard spectral state with Chandra X-ray Observatory observations on 2017 February 6 (MJD 57790, Miller et al. 2017) and International Gamma-Ray Astrophysics Laboratory (INTEGRAL) observations on 2017 February 10 (MJD 57794, Del . The source was then seen transitioning to the hard-intermediate state for some time with Neil Gehrels Swift Observatory (Swift) observations on 2017 March 27 and April 2 (MJD 57839 and MJD 57845, Armas Padilla & Munoz-Darias 2017), and then returning to the hard state after a failed-transition outburst to the soft state on 2017 May 5 and 11 (MJD 57878 and MJD 57884, Bassi et al. 2017), as was also the case in the 1993 event (Revnivtsev et al. 1998). Bassi et al. (2019) studied the HID of the source and found that it had three softening events when the source transitioned from the hard to the hard-intermediate state. Along with the three softest points (MJD 57854.2, 57895.9 and 57960.7), we consider all the dates with hardness ratio 0.7, which lies in the range of 2017 July 6 and August 13 (MJD 57940-57978) as the hard-intermediate state.
The source was found to be one of the 'outlier' BHXBs (Bassi et al. 2019) in the radio/X-ray correlation plane (which are radio fainter by 1-2 orders of magnitude, and tend to have a steeper correlation index, with β ∼ 1.4, e.g., Corbel et al. 2004;Coriat et al. 2011;Gallo, Miller & Fender 2012). A type-C quasi-periodic oscillation (QPO) was also detected in the hard state (Bharali et al. 2019), and signatures of a hot and dense accretion disk wind (with terminal velocity ∼2000 km s −1 ) were observed (Cuneo et al. 2020). From the broadband spectral fitting of the source, the irradiated accretion disk was found to dominate the optical emission, while a hint of an excess near-IR emission above the prediction of the irradiated disk model was observed, likely due to synchrotron emission originated in the jet (Rout et al. 2021).

System parameters of GRS 1716−249
The system parameters of GRS 1716−249 are not well constrained. From the 1993 outburst, della Valle et al. (1994) proposed that the system contains a low-mass main-sequence star with spectral type K (or later), at a possible distance between 2.2 kpc (lower limit obtained from the equivalent width of the NaD absorption lines) and 2.8 kpc (upper limit based on an incorrect maximum luminosity of an X-ray transient). But in light of several arguments we explore in Section 4.2.2., we find that the estimated upper limit of 2.8 kpc is not a reliable constraint for its distance. Masetti et al. (1996) discovered superhumps in the lightcurve (although these could also be due to irradiation modulation, see Section 4.2.1 for a discussion). Assuming that the donor is a main sequence star, they estimated the companion star mass to be ∼1.6 M and inferred an orbital period of ∼0.6127 days or ∼14.7 hrs for a Roche lobe filling star. Then they used the maximum mass ratio criterion for having superhumps, which is about 3:1, and proposed that the mass of the accreting compact object is > 4.9 M , hence classifying it as a black hole. They also suggested that a 1.6 M main sequence star at 2.4±0.4 kpc would exceed the quiescent luminosity of the binary substantially (although it is important to note that the quiescent luminosity limit of the source was not confidently known, see Section 3.6). Despite all the crude assumptions employed, these limits on the mass of the compact object and the distance to the source have been used for all subsequent studies on the source, until this paper.
During the 2016-2017 outburst, GRS 1716−249 was extensively studied in the X-ray wavelengths. Tao et al. (2019) used spectral fits of three NuSTAR and Swift datasets in its hard-intermediate state, and constrained the black hole mass to be < 8 M at a 90% confidence level under the assumption that the distance to the source is 2.4±0.4 kpc. Using the same assumption, they also inferred the inclination angle of the inner disk  (MJD 57874.7) with the 2-m LCO telescope in the i -band with 200s exposure time. Previously, a lower-resolution optical finding chart during outburst is available in the V band (Masetti et al. 1996). The target is indicated with hash mark in both the panels. The right panel shows the quiescent optical finding chart (MJD 58388.4) with the 2-m LCO telescope in the i -band with 100s exposure time (an image taken under excellent conditions, with seeing of 0.82 arcsec). The counterpart is just 1.6 arcsec from a star of similar magnitude to the north of GRS 1716-249. to be in the range of 40 • -50 • by performing joint modelling of the continuum and the reflection components. An analysis of the broadband (1-78 keV) X-ray spectra of the source taken by NuSTAR and Swift constrained the accretion disk density parameter of GRS 1716−249 to be in the range of 10 19 -10 20 cm −3 (Jiang et al. 2017). Recently, the black hole mass was claimed to lie in the range of 4.5-5.9 M according to a two-component advective flow (TCAF) model (Chatterjee et al. 2021), although this method uses model-dependent spectral fitting of the source to obtain these values.
In this paper, we present a detailed multi-wavelength study of GRS 1716−249 during its 2016-2017 outburst, with particular focus on its UV/optical/IR emission to investigate the physical mechanisms contributing to the emission in these wavebands, and reveal the system parameters of the source. In Section 2, we describe in detail the observations and the analyses of the data used for this study. In Section 3, we present the characteristics of the outburst using various tools like the light curves, variability of the source during the peak of the outburst using fractional rms values, the optical/UV spectra of the source, the broadband spectral energy distributions (SEDs), the colour-magnitude diagrams to study the colour evolution of the source during the outburst, and the optical/UV/X-ray correlations to explore the various emission mechanisms. We also report longterm (∼ 10 years) monitoring of the source and discuss its quiescent optical magnitude, which is important as the optical brightness of BHXBs in quiescence has minimal contribution from the accretion disk and is dominated by the companion star (Chevalie et al. 1989). In Section 4, we interpret and discuss our results, including the implications of our analyses on the system parameters of the source, and present new estimates for the distance to the source. Finally, we present our conclusions in Section 5. We monitored GRS 1716−249 during its 2016-2017 outburst extensively with the Las Cumbres Observatory (LCO) between 2017 January 28 and October 21 (MJD 57781-58046). Observations were made using the 1 m LCO telescopes at Siding Spring Observatory, Australia, Cerro Tololo Inter-American Observatory, Chile, and the South African Astronomical Observatory (SAAO), South Africa, as well as the 2 m Faulkes Telescopes at Haleakala Observatory, Maui, Hawai'i, USA and Siding Spring Observatory, Australia. The source was also monitored during quiescence, before and after the 2016-2017 outburst, for 11 years since 2006 February, as part of an on-going monitoring campaign of ∼50 low-mass Xray binaries (LMXB) coordinated by the Faulkes Telescope Project (Lewis et al. 2008;Lewis 2018).
Imaging data were primarily taken in the Sloan Digital Sky Survey (SDSS) g , r , i and the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) Y -band filters, with some data also taken in Bessel B and V -bands. The data were initially processed using the LCO BANZAI pipeline (McCully et al. 2018). The multi-aperture photometry on the reduced data was performed using "X-ray Binary New Early Warning System (XB-NEWS)", a real-time data analysis pipeline that aims to detect and announce new X-ray binary outbursts within a day of the first optical detection of an outburst (e.g. Pirbhoy et al. 2020;Goodwin et al. 2020). The XB-NEWS pipeline downloads new images of all targets of interest from the LCO archive along with their associated calibration data and performs several quality control steps to ensure that only good quality images are analysed. XB-NEWS then computes an astrometric solution for each image using Gaia DR2 positions 1 , performs aperture photometry of all the stars in the image, solves for zero-point calibrations between epochs (Bramich & Freudling 2012), and flux calibrates the photometry using the ATLAS All-Sky Stellar Reference Catalog (ATLAS-REFCAT2, Tonry et al. 2018). The pipeline also performs multi-aperture photometry (azimuthally-averaged PSF profile fitting photometry, Stetson 1990). Light curves are produced in near realtime. If the location of the source is well-known, but the source is fainter than the formal detection threshold, the pipeline performs forced photometry on the position. Magnitude errors larger than ∼ 0.25 mag are considered as marginal detections and are not included in our study.
We detected the source during outburst in a total of 192 images between 2017 January 28 (MJD 57781) and 2017 October 21 (MJD 58046), generally at a cadence of every 2-3 days during the brighter phase of the outburst, and every ∼75 s for the high cadence images taken on 2017 May 9 (MJD 57882). A detailed observation log containing information about the LCO epochs, filters and magnitudes, is summarized in the Appendix (see Table A1). From our XB-NEWS optical analysis, the accurate optical position of the source was found to be RA: 17:19:36.917 and Dec: -25:01:04.20 (J2000), consistent with the VLBI coordinates to within 0.1 arcsec (Atri et al. 2019). The optical finding charts in the iband during both outburst and quiescence, are shown in Fig. 1. The systematic error in the position measurement is small ( 0.3 arcsec) and has better precision than previously reported optical measurements (∼1 arcsec in della Valle et al. 1994). To convert the multi-aperture photometry magnitudes obtained from the XB-NEWS pipeline to the intrinsic de-reddened flux densities, we use the absorption column density N H = (0.70±0.01) ×10 22 cm −2 , as reported in Bassi et al. (2020). Using the relation between optical extinction and hydrogen column density (Foight et al. 2016), the V -band absorption coefficient is inferred as A V = 2.44 ± 0.11 mag. There are different determinations of the relation between ex-1 https://www.cosmos.esa.int/web/gaia/dr2 tinction and hydrogen column density in the literature (see for e.g. Guver et al. 2009;Watson 2011;Willingale et al. 2013), but we use the Foight et al. (2016) value as they provide the most recent estimates using updated abundances and hence are likely more reliable. We note that Bahramian et al. (2015) also arrive at a similar relation using the updated abundances, while including X-ray binaries in their sample. This leads to a color excess of E(B − V )∼0.8 mag (assuming a mean value of A V /E(B − V )∼3.1 for the diffuse interstellar medium, Fitzpatrick 1999), which is consistent with the historical value of E(B − V )∼0.9±0.2 mag (della Valle et al. 1994) obtained based on multiple lines of reasoning. The wavelength-dependent extinction terms, used for de-reddening in other bands, are obtained from the extinction curve of Cardelli et al. (1989).

Archival optical data
We also use the archival data of the source obtained in the G spectral filter with the Gaia telescope 2 during the recent outburst. Gaia first detected the source on 2017 January 27 (MJD 57780.8) at G=16.44. Prior to that, the last observation it had on 29 October 2016 (MJD 57690) was a non-detection (typical detection limit of Gaia is ∼ 20.7 mags, Brown et al. 2016). Gaia detected GRS 1716−249 on 13 days during the outburst, with the last detection on 2017 September 23 (MJD 58019). We use these public data in Fig. 2, while studying the optical light curve of the source.

Archival historical optical data
To compare the 2016-2017 outburst of the source with its discovery outburst from 1993, we include the simultaneous optical detections taken on 1993 October 8 (MJD 49268) in the B, V and R filters as 17.7±0.1, 16.7±0.1 and 16.0±0.1 mags, respectively (della Valle et al. 1994). We also use the historical B and V -band observations from Masetti et al. (1996). We use these data in the spectral energy distribution study (see Section 3.3) and the detailed analysis of the evolution of the source through the colour-magnitude diagram (see Section 3.5).

VISIR mid-IR observations
We acquired targeted observations of GRS 1716−249 with the Very Large Telescope (VLT) in mid-IR wavelengths on three nights during the 2016-2017 outburst, using the VLT Imager and Spectrometer for the mid-IR  Table 1. The data were reduced using the VISIR pipeline in the gasgano environment. We combined the raw images from the chop/nod cycle and performed aperture photometry in IRAF using a large enough aperture to minimize the effect of small seeing variations on the fraction of flux in the aperture (the method is the same as that used in Baglio et al. 2018). To flux calibrate the photometry and estimate the flux density of the source, we used all the standard star observations taken within one month of the observation night in the same filter during clear sky conditions. All the standard stars used are listed in the final column of Table 1. At mid-IR wavelengths, the zero point corrections rarely vary much. In fact, we found that the ADU/flux conversion factor measured from different standard star observations within a month varied only by 5%-10% when the airmass is less than 1.5. For all the filters for which we had only one standard star available within one month of the observation, we use an error of 5% in the ADU to calculate the uncertainty on the flux density of the source.
The source was only detected on one date, 21 April, in the J8.9 filter, with a magnitude of 18.24±0.19, or a flux density of 3.22±0.59 mJy (the detection has a signalto-noise ratio of 6.2). Although the photometric error for the detection is small, the uncertainty in the flux density is increased due to systematic errors arising from the limited number of available standard stars within a month of the detection. On the other two dates when the source was not detected, we derive 3σ upper limits from the root mean square (rms) in a region centred on the position of GRS 1716−249. The closest WISE catalogue star is 12 arcseconds away from the position of GRS 1716−249 (outside the field-of-view of VISIR), with a flux density of 2.24 mJy at 12µm.

REM near-IR observations
We observed GRS 1716−249 in the near-IR wavelengths (J, H and K bands, one filter at a time) with the REMIR camera mounted on the Rapid Eye Mount (REM; La Silla, Chile) telescope between 2017 February 8 and October 1 (MJD 57792-58027). For each epoch, the reduction of the images was performed by subtracting the sky contribution; this was obtained as the median of 5 misaligned exposures of 60s and 30s in the J and H filter, respectively, and of 10×15s exposures in Kband. Once the sky was subtracted, we registered and averaged the exposures to enhance the signal to noise.
We performed aperture photometry on each reduced image using IRAF. The magnitudes were then calibrated against a group of five 2MASS reference stars in the field. We note that a 2MASS star is observed at a distance of ∼ 3.5 from the target in the REM images. Considering the spatial resolution of REMIR (1.22 /pixel) and the seeing, the two stars are therefore blended together in all images. Under the reasonable hypothesis that the 2MASS star is not variable, we subtract the contribution of the 2MASS star from the flux extracted with our analysis to build spectral energy distributions. The magnitudes of the 2MASS star are tabulated in the 2MASS catalogue (J = 15.33 ± 0.06; H = 14.46 ± 0.06; K = 14.16±0.08). To double check, we found and downloaded archival J-band images of the field taken with the SOFI instrument at the New Technology Telescope (NTT; La Silla, Cile) during quiescence in 1999 (July 5 and 7; Program ID: 63.H-0232). The J-band magnitude of the 2MASS star, after calibration, is J = 15.38±0.05, entirely consistent with the value reported in the 2MASS catalogue (which suggests that source has probably been stable over the years). We tabulate the REM epochs, filters and magnitudes in the Appendix (see Table A2).

Archival near-IR observations
We use the archival near-IR photometric observations of GRS 1716−249 during the outburst. Bassi et al. (2020) reported near-IR detections of the source obtained with the Rapid Eye Mount telescope (REM) on 2017 February 9 (MJD 57793) of J=14.18±0.22, H=13.81±0.14 and K = 13.84±0.29 and 13.59±0.16, with exposure times of 300s, 150s, 75s and 75s, respectively (the same observation is also included in the dataset presented in Section 2.2.2). Later, Joshi et al. (2017) observed the source with the Mount Abu 1.2 meter telescope and the Physical Research Laboratory (PRL) near-IR Imager/Spectrograph, and reported near-IR magnitudes on 2017 March 20 (MJD 57832) of J = 14.3, H = 14.0 and K S = 13.7, with typical errors of 0.1 magnitude, for a total integration time of 15, 15 and 17.5 minutes, respectively. These magnitudes are consistent with those derived in our analysis of the REM data during the 2016-2017 outburst, before the subtraction of the contribution from the nearby 2MASS star (see Section 2.2.2 for a detailed discussion).
We also use the historical IR data of the source from its mini-outburst in 1994, taken on 1994 July 8 (MJD 49541) in the J and K filters of 16.2±0.3 and 15.5±0.3, respectively (Chaty et al. 2002), especially for the spectral energy distribution study (see Section 3.3).

Archival Swift/UVOT Observations
We gathered publicly available Swift UV/Optical Telescope (UVOT) pointing observations of the source during its entire outburst from the NASA/HEASARC data center. We used the pipeline processed images and obtained the magnitude of the source using the uvotsource HEASOFT routine, with an aperture of 5 arcsec centered on the source. An empty region close to GRS 1716−249 was chosen as the background region. We select only those 81 pointings where the source flux estimate is at least 5σ above the sky background. Although the UVOT observations of this source were available in all the filters, most of the significant and usable detections were found to be in the V (0.546 µm), B (0.439 µm) and U (0.346 µm) bands, with a smaller amount of detections in the U V W 1 (0.260 µm) bands.
Similar to the optical flux values, we de-reddened the UV data. We use the absorption column density N H = (0.70±0.01) ×10 22 cm −2 , reported in Bassi et al. (2020), the generic relation between optical extinction and hydrogen column density (Foight et al. 2016), and the wavelength-dependent extinction terms using the extinction curve of Mathis (1990).

Archival radio detections
We use the radio observations of the source during its outburst with the Karl G. Jansky Very Large Array (VLA; 5.25, 7.45, 8.8 and 11.0 GHz), Australia Telescope Compact Array (ATCA; 5.5 and 9.0 GHz) and Australian Long Baseline Array (LBA; 8.4 GHz) as reported in Bassi et al. (2019) and Atri et al. (2019). Radio detections of the source are available for 2017 February 9 and 11 (MJD 57993 and 57795), April 22 (MJD 57865), August 12 and 13 (MJD 57977 and 57978). We use all the radio detections for which more than one quasi-simultaneous (within 24 hours) optical/UV measurement is available, to study the broadband spectral energy distribution (SED) of the source (see Section 3.3).

Archival data from X-ray telescopes
We acquired X-ray monitoring data of GRS 1716−249 from the Swift/BAT and Swift/XRT telescopes. Swift/BAT has observed the source almost daily from 2016 December 1 (MJD 57723) in the 15−50 keV flux range. We extracted the daily average light curve data of this source from the Swift/BAT transient monitor 3 (Krimm et al. 2013). To convert the count-rates to flux, we used the hydrogen column density N H = (0.70±0.01) ×10 22 cm −2 and a photon index of Γ = 1.68±0.01, as reported in (Bassi et al. 2020). Swift/XRT observations were made every few days between 2017 January 28 (MJD 57781) and 2017 October 20 (MJD 58046), mostly in the window timing mode, with target IDs 34924 and 88233 (see Table 1, of Bassi et al. 2019, for a detailed observation log). We used the on-line Swift/XRT products generator 4 (Evans et al. 2007(Evans et al. , 2009) to extract the 2−10 keV count rate of GRS 1716−249 from each XRT observation, after correcting for instrumental artefacts.
GRS 1716−249 was also observed with Nuclear Spectroscopic Telescope Array (NuSTAR) during the 2016-2017 outburst. We calculate the NuSTAR flux density of the source using NUPIPELINE V0.4.6 in HEASOFT V6.25, with the calibration file version v20171002. Both Science Mode and Spacecraft Mode data were considered. The flux was calculated in the Swift BAT energy band (15-50 keV) using the best-fit spectral models provided in Jiang et al. (2017). A relativistic disk reflection model with a variable disk density parameter was used (Ross & Fabian 2007). Detailed descriptions of spectral modelling can be found in Jiang et al. (2017).
We also gathered daily X-ray monitoring data of GRS 1716−249 from MAXI/GSC 5 (Matsuoka et al. 2009) in the 2−20 keV range covering the complete outburst.

Multi-wavelength light curve
The light curves of the entire outburst are plotted in Fig. 2 including data in near-IR (REM), optical (LCO and Gaia), UV (UVOT) and X-ray (NuSTAR, Swift XRT and BAT, MAXI) wavelengths.
The first optical detection of GRS 1716−249 during the 2016-2017 outburst was obtained by LCO on MJD 57781, when the optical magnitude was already bright, with i = 15.97 ± 0.01. Since that time we regularly monitored the source in i , g , r and Y bands until its optical emission faded back to quiescence. There are also some scattered observations taken in the B and V bands. A zoom-in of the optical light curve during the peak of the outburst between 2017 January 28 and May 27 (MJD 57800-57900) is shown in Fig. 3a. The general trend of the LCO light curve is an almost constant plateau in all the optical bands during the outburst, until ∼ 2017 May 27 (MJD 57900), followed by a slow and steady decline to quiescence until 2017 October 20 (MJD 58046). The same behaviour was also observed in the Gaia optical light curve in G-band. At the end phase of the decline, we find a small amplitude optical brightening of the source, with i magnitudes changing from 17.80 ± 0.01 on 2017 October 7 (MJD 58033), to 17.45 ± 0.01 on October 15 (MJD 58041), and then again back to 18.10 ± 0.01 on October 18 (MJD 58044).
The optical/UV light curves obtained with Swift/UVOT in the U , B, V , U V W 1, U V W 2 and U V M 2 bands, show a similar outburst profile as the 5 http://maxi.riken.jp/top/lc bh.html LCO and Gaia light curves. The complete outburst, including the decay towards quiescence, is well-covered by the U -band data. The brightening observed in optical wavelengths during the decline of the outburst is not evident in the UV data.
The near-IR REM light curve is approximately constant, with some scatter and flickering, until 2017 June 21 (MJD 57925). A zoom-in of the near-IR light curve during the peak of the outburst is shown in Fig. 3b. After MJD 57925, the flux starts to show a decreasing trend in all bands until ∼ 2017 September 4 (MJD 58000), when the flux experiences a plateau that lasts until the end of the observations. This behavior is similar to that of the higher energy light curves.
The X-ray light curves in both the hard (Swift/BAT and NuSTAR) and the soft (Swift/XRT and MAXI) energy ranges follow the morphology of the near-IR/optical light curve, indicating a correlated behaviour, which is explored in detail in Section 3.4.  Fig. 3a shows the zoom-in of the optical light curve during the peak of the outburst. On longer timescales, the four LCO optical filters (g , r , i and Y band) are clearly correlated. We also took higher cadence optical observations on 2017 May 9 (MJD 57882; 15 detections in ∼17.5 mins with a time resolution of ∼ 75 seconds) of the source with LCO i -band. The optical fractional rms deviation in the flux on such short timescales (minutes; i.e. a frequency range of 0.0010-0.013 Hz) during the hard state, evaluated following the method described in Vaughan et al. (2003) and Gandhi et al. (2010), is found to be 1.3±0.4%, reflecting on a very weak variability.

Optical variability
The observed rms is substantially lower than the optical fractional rms of BHXBs like GX 339-4 and Swift J1357.2-0933 in the hard state and V404 Cyg in the flaring state, which are ∼ 5-20% on similar and shorter timescales (Gandhi 2009;Gandhi et al. 2010;Cadolle Bel et al. 2011;Gandhi et al. 2016;Paice et al. 2019). The variability seen in GRS 1716−249 is similar to the lower optical fractional rms values of ∼3-5% seen in the hard accretion states of Swift J1753-0127, XTE J1118+480 and MAXI J1535-571 (Gandhi 2009;Hynes et al. 2009;Baglio et al. 2018). Such variability is also observed in the soft accretion states of GX 339-4 and GRO J1655-40 (Hynes et al. 1998;O'Brien et al. 2002;Cadolle Bel et al. 2011), when the accretion disk dominates the emission. Fig. 3b shows the zoom-in of the near-IR light curve during the peak, where the source is observed to be varying by ∼1 magnitude. We also observe a possible small amplitude (∼0.5-0.6 mag) flare happening in J and K bands between MJD 57887 and 57891. However, no corresponding activity is observed in H-band, and the lack of time-resolved data during these days makes it difficult to study this event further. The fractional rms deviation in the infrared flux of GRS 1716−249 during the peak of the outburst on longer timescales (days/weeks; a frequency range of 5.8E-6 -8.7E-8 Hz), after removing the contribution from the blended star, is measured to be 20.69±2.34%, 10.92±4.86% and 34.43±4.24% in the J, H, and K bands, respectively. Hence GRS 1716−249 is variable in the near-IR band. Although the coverage of the outburst in the mid-IR range is scarce, the detections and the upper limits (see also Table 1 and Section 2.2.1) also point to a variable mid-IR emission, with the flux density spanning from < 1.4 mJy to 3.2 ± 0.6 mJy at 8.7 µm.

Infrared variability
A similar increase of fractional rms deviation in the flux towards longer wavelengths in the optical/IR wavelength range, is also seen in other BHXBs. For example, the rms is often 10-20% or higher in the near-IR regime in the hard accretion states, as seen in sources like XTE J1550-564, GX 339-4 and MAXI J1820+070 (Curran & Chaty 2013;Vincentelli et al. 2018;Tetarenko et al. 2021). In the mid-IR regime, the fractional rms increases further, with for example rms ∼15-22% in MAXI J1535-571 at a similar time resolution, which supports a jet origin to the variability on these timescales (Baglio et al. 2018). In XTE J1118+480, the spectrum of the rms variability is consistent with a power law of spectral index α = −0.6 from optically thin synchrotron radiation, spanning near-IR to X-ray (Hynes et al. 2003. In the hard accretion state of MAXI J1820+070, the fractional rms (in a larger integrated frequency range; 10 −4 -50 Hz) decreases monotonically with increasing wavelength, from tens of per cent in the optical/near-IR, to 2-8% at radio frequencies (Tetarenko et al. 2021). Other timing properties such as the frequency of the break in the power spectrum was also seen to vary smoothly with wavelength from optical to radio, with time lags between bands increasing at longer wavelengths. One interpretation is that although the fractional variability increases from optical to IR due to an increase in the jet contribution, the fractional rms drops again as it approaches radio wavelengths, because the variability in the jet dominated bands gets more smoothed out by the larger size scale of the emitting region at the longer wavelengths.

Spectral energy distribution
We construct the optical/UV spectra as well as the broadband SEDs of GRS 1716−249 in the hard (Fig. 4a) and hard-intermediate states (Fig. 4b), to illustrate the peculiar multi-wavelength characteristics of the source. In Fig. 4 we plot the optical/UV spectra of GRS 1716−249 in both the hard (Fig. 4a) and hard-intermediate states (Fig. 4b) of the outburst. We use quasi-simultaneous observations obtained within 24 hours, and convert the magnitudes to de-reddened fluxes as described in Section 2.1 for optical LCO magnitudes and Section 2.3 for archival Swift/UVOT observations. In both the hard-intermediate and hard states, the SEDs are found to be smooth up to the U V W 1-band, with a shallow peak around the g band.
We use the available information to constrain the intrinsic optical/UV spectral index by fitting the function S ν ∝ ν α , where S ν is the flux density, ν is the frequency and α is the spectral index. We obtain an average U -i spectral index of α U −i = −0.1±0.3 across The de-reddened radio/mid-IR/near-IR/UV spectrum when quasi-simultaneous (within 24 hr) data were available. The mid-IR to radio spectral index measured from the VISIR mid-IR detection on MJD 57864.4 and LBA radio detection at 8.4 GHz on MJD 578865.7 is found to be α=0.13±0.03. The mid-IR upper-limit on MJD 57837 in the J8.9 filter is plotted as an inverted triangle, to show the mid-IR variability of the source. We also plot for reference the historical optical (MJD 49268, della Valle et al. 1994) and near-IR (MJD 49541, Chaty et al. 2002) SEDs from its discovery outburst in 1993/1994 (in grey dotted lines). the spectra. Generally, a negative slope ∼ −0.7 (e.g. Gandhi et al. 2011) is expected if there is a jet present with an optically thin synchrotron spectrum dominating the near-IR/optical regime. This value can be even more negative, as seen in cases like Swift J1357.2−0933 where the quiescent optical/mid-IR SED has a powerlaw index of −1.4, arising from a weak jet (Shahbaz et al. 2013). Although optically thick, self-absorbed synchrotron emission from a jet can produce slope ∼ −0.1, there are only a few cases in which such emission extends to higher frequencies like optical (e.g. Russell et al. 2013;Maitra et al. 2017). A positive slope SED is expected (with spectral index ∼ 1) if the optical emission is dominated by the blackbody from the outer accretion disk (e.g. Hynes 2005). For a viscously heated disk, α ∼ 0.3 is expected, turning to a steeper slope 0.3 < α < 2.0 at lower frequencies (e.g. Frank et al. 2002). Very often, a combination of all the processes can result in an intermediate slope.
For comparison, during the hard state, the spectrum constructed for GRS 1716−249 on 2017 February 22 (MJD 57806), has an α Y −g ∼ 0.5, while during the hard-intermediate state on 2017 April 11 (MJD 57854), we find α Y −B ∼ 0.7 (See Fig. 4a and 4b). This suggests the optical spectra are probably dominated by an accretion disk, but it is unlikely for the UV/optical emission to solely originate from the Rayleigh-Jeans part of a single-temperature blackbody spectrum. The quasi-simultaneous, de-reddened NIR/optical/UV spectrum of GRS 1716-249 on MJD 57864 (blue circles). Superimposed is the fit of the spectra with a single-temperature black body (black line). We also overplot the quasi-simultaneous mid-IR detection of the source for comparison (red diamond). We show that the NIR/optical/UV part of the spectrum is qualitatively well represented by the single temperature black body, while the mid-IR emission is comparatively brighter probably due to an additional contribution from the jet.
To investigate the issue further and disentangle the emission processes in the optical and UV regime, we include the available IR and radio data and construct the broadband spectrum (see Fig. 4c) and the SED (see Fig. 6) of GRS 1716−249 with quasi-simultaneous (within 24 hours) optical (with LCO), near-IR (with the Mount Abu 1.2 meter telescope and REM), mid-IR (with VISIR) and radio (with ATCA, LBA and VLA) data. The broadband spectrum has a positive slope in the near-IR regime, which flattens in the optical, and gets fainter in the UV wavelengths. We fit the spectrum on MJD 57864 with a single temperature black body curve (See Fig. 5), and find that the NIR/optical/UV part of the spectrum is qualitatively well represented by a black body model with a temperature ∼10500±200 K, while the mid-IR emission is comparatively brighter. This suggests that the overall shape of the IR/optical/UV spectra together are as expected for the outer accretion disc. The steeper slope in the NIR is the Rayleigh-Jeans limit of the blackbody with the lowest temperature. The fainter UV emission suggests that the viscous disk does not dominate in these wavelengths (as the UV emission does not keep rising with alpha ∼0.3). Instead, the irradiated disc most likely dominates the emission, as a peak is seen around g'-band, with UV flux densities being slightly fainter. This is similar to seen in other BHXBs where the irra-diation bump peaks in the optical, with the UV slightly fainter (e.g. Hynes 2005). The historical optical (taken in 1993) and near-IR (taken in 1994) data from the discovery outburst, plotted in Fig. 4c for a comparison with the current outburst, show a similar brightness profile in optical, although not as flat as the recent data. In the near-IR wavelengths, although it was fainter during the historical outburst, it shows a similar shape with a positive slope.
The single mid-IR detection of the source, obtained on 2017 April 21 (MJD 57864), is significantly brighter than what is expected from the disk alone (see Fig. 5 and Fig. 6), and can probably be attributed to synchrotron emission from a compact jet during the outburst. The mid-IR to radio spectral index measured from the VISIR mid-IR detection on MJD 57864.4 and the LBA radio detection at 8.4 GHz on MJD 57865.7 (a separation of 1.3 d) is found to be α=0.13±0.03. This is slightly more positive than (but consistent with to within < 2σ), the reported radio spectral indices of α=-0.15±0.08 and α=-0.07±0.19, at the beginning (2017 February 9, MJD 57793.8) and close to the end of the outburst (2017 August 12, MJD 57977.3), respectively (Bassi et al. 2019). Therefore, the radio to mid-IR spectrum is consistent with a flat or slightly inverted spectrum coming from a compact jet. We also estimate the spectral index between the mid-IR detection and quasi-simultaneous Xray spectrum as α = −0.26. If the mid-IR emission arises from optically thin synchrotron, then its extrapolation to X-ray is much fainter than the observed X-ray power law index, which implies that the synchrotron jet does not contribute much to the X-ray flux.

Multi-wavelength correlations
Another tool for disentangling the emission processes in BHXBs during outburst is multi-wavelength correlations. We study the quasi-simultaneous multiwavelength correlation of GRS 1716−249, using dereddened optical and UV fluxes as a function of the soft X-ray (2-10 keV) fluxes from Swift/XRT, and hard Xray (15-50 keV) fluxes from Swift/BAT and NuSTAR; whenever the X-ray fluxes are obtained within a day of the optical or UV observations.

Optical versus X-ray correlations
For the optical versus soft X-ray (2-10 keV) correlation study, we use Swift/XRT flux in the 2-10 keV range for X-rays, and de-reddened optical i -band flux density (as i -band had the best coverage amongst all optical filters). We choose all the points for which we have quasi-simultaneous data (i.e., data obtained within 24 hours; see Fig. 7a). While the hard state values follow one single correlation, the hard-intermediate state values (shown in the plot as coloured points with noncircular symbols, where the MJD values are indicated) are generally seen to lie on the lower side of the correlation. This is in agreement with previous studies where comparatively less optical emission is observed during the state transition and soft state (eg. Jain et al. 2001;Corbel & Fender 2002;Russell et al. 2006;Coriat et al. 2009). Generally, this is thought to be due to a weak jet component to the optical emission, which usually fades during the transition from the hard to hard-intermediate state (eg. Cadolle Bel et al. 2011;Baglio et al. 2018) and recovers when a BHXB returns to the hard state (e.g. Corbel et al. 2013;Kalemci et al. 2013;Russell et al. 2013). Another reason could be a weak disk-blackbody component which can contribute towards the X-ray luminosity in the 2-10 keV energy range during the hard-intermediate state (Capitanio et al. 2009;Alabarta et al. 2020).
The correlation is found to be significant (Pearson correlation coefficient = 0.84, p-value = 2.5 × 10 −10 ). The best-fit slope for the correlation in hard state using the orthogonal distance regression method of least squares is 0.41±0.03. The observed slope suggests an X-ray irradiated accretion disk (van Paradijs & McClintock 1994), with possibly some contribution from the viscous disk (Russell et al. 2006). But we note that the scaling relation of van Paradijs & McClintock (1994) depends on an assumed geometrical configuration and is not as simple and straightforward. Recent studies have showed that the slope of the correlation can differ depending on the origin of the emission at different regimes (e.g. Coriat et al. 2009;Tetarenko et al. 2020), and irradiation from a hot dense accretion disk wind may also cause a slight distortion of the scaling relation (Cuneo et al. 2020, see also Section 4.1 for a detailed discussion).
In addition, we note that there is a slight hint of the correlation flattening at the fainter end of the luminosity ranges. As an alternative explanation, we attempted to fit the data with a broken power law (keeping the break luminosity as a free parameter). The correlation obtained were found to have a steeper slope (∼1.1) at the brighter end, and a shallower slope (∼0.2) at the fainter end, implying that the viscous disk could play a role at the lower luminosities (the break luminosity was found to be ∼ 2.5 × 10 −9 erg s −1 cm −2 ). A prominent role of the viscous disk in the fainter part of the outburst is also hinted at by our color evolution analysis (see Section 3.5).
The hard X-ray emission (15-50 keV) of the source follows a power-law correlation with the optical flux.
To check this correlation, we used hard X-ray data in the 15-50 keV range from Swift/BAT telescope and NuSTAR with de-reddened i -band flux density obtained from LCO (see Fig. 7b). The correlation with a power-law index of ∼0.54±0.04 is found to be significant (Pearson correlation coefficient = 0.93, p-value = 1.2 × 10 −13 ). The difference between hard state and hard-intermediate state here is subtle, as probably the hard X-ray flux is also fading slightly in this state compared to the hard state, such that the optical and the hard X-ray flux are both fainter, maintaining the correlation.

U -band versus X-ray correlations
The U -band/soft X-ray correlation is plotted in Fig. 7c with U -band detections from the UVOT telescope, and simultaneous soft X-ray data from Swift/XRT in the 2-10 keV energy range. The correlation is significant (Pearson correlation coefficient = 0.94, p-value = 2.6 × 10 −12 ) and the power-law index of the UV/X-ray correlation is found to be ∼0.49±0.03, which is consistent with an irradiated accretion disk (van Paradijs & McClintock 1994). Similar to the optical/X-ray correlation, the hard-intermediate state values were seen to have weaker U -band emission in comparison to the hard state. A correlation was also observed between the U -band and the hard X-ray emission (15-50 kev, from Swift/BAT telescope and NuSTAR, see Fig. 7d), with high significance (Pearson correlation coefficient = 0.95, p-value = 1.2 × 10 −15 ), and a similar slope of ∼0.51±0.03. There is a hint of the correlation flattening at the lower luminosity end, as also seen in the case of the optical/X-ray correlations. But due to the lack of fainter data points, and the large uncertainties associated with it, a conclusive result regarding a broken power law can not be obtained. But we note that a shallower correlation could arise due to the emergence of a viscous disk at the end of the outburst (see also Section 3.5).

Colour-magnitude diagram
The colour-magnitude diagrams (CMDs) are plotted in Fig. 8, using four different filter combinations using the i , g , U , B and V filters. We adopt the sin- Figure 8. Colour-magnitude diagrams, adopting four different filter combinations. For each combination, the bluer filter is on the y-axis. Bluer colours (greater spectral index) are to the left, redder colours (lower spectral index) are to the right. A simple model of a single temperature blackbody heating up and cooling, used to approximate emission from reprocessing on the disk, is denoted by the black line labelled with the temperature values, in each panel (see text). The U − g CMD (a) is the combination showing the shortest wavelengths; the reprocessing model is a poor approximation of the data, the viscously heated disk likely plays a strong role at the fainter epochs. The B − V CMD (b) also includes some data from the 1993 outburst (della Valle et al. 1994;Masetti et al. 1996); the reprocessing model approximates most of the data well (brightest epochs in both outbursts). The g − i CMD (c) is the combination with the longest wavelengths; the reprocessing model is close to the data at bright epochs but not during the outburst fade, when the viscous disk likely dominates. The U − i CMD (d) shows the widest wavelength range; again reprocessing can describe most of the brightest epochs, not the fainter epochs. gle temperature blackbody model of Maitra & Bailyn (2008), described in detail in  to study the colour evolution of X-ray binaries during outbursts, which approximates the emission from the X-ray irradiated outer accretion disk. A colour change is expected due to the evolving temperature of the irradiated disk, which is assumed to have a constant emitting surface area. While at high temperatures the optical emission is expected to originate in the Rayleigh-Jeans tail of the blackbody, at lower temperatures it originates near the peak of the blackbody curve. The blackbody temperature of the model depends on the intrinsic colour and the interstellar extinction.
The normalization of the model depends on the accretion disk radius (estimated using the known orbital period, mass of the companion star and the mass of the black hole from the literature), the distance to the source, inclination angle, the disk filling factor, disk warping and the fraction of disk that is reprocessing the X-rays. As many of these parameters are not certain, we choose a value of the normalization that best describes the trend in the data. In particular, we fix the normalization using the B-V CMD (see Fig. 8b) as this filter combination has the most amount of data, and is less affected by the jet emission (if present), being from the bluer wavelengths. We find that the data do not completely agree with the single temperature black-body model, which indicates that more than one component is likely to be present. The disk temperature was roughly seen increasing from ∼7,000 K to ∼12,000 K, as expected during outbursts when hydrogen in the disk is typically ionized.
We find that the data in the B-V CMD, which includes some data from the 1993 outburst (data from della Valle et al. 1994;Masetti et al. 1996), generally follows the expected trend between colour and magnitude, with scatter of ±0.1 mag in colour. These data were from the brightest epochs in both outbursts. We adopt the same normalization that this provides to the other three filter combinations. We note that there is an uncertainty on the normalization due to the scatter in the colour, but it should not be larger than ∼0.1 mag in colour. Assuming that the same normalization can be applied to the other filter combinations, we can investigate deviations from the blackbody model as a function of wavelength combination. For these filter combinations, we find that the brightest epochs show data close to the blackbody model, but at lower luminosities there are significant deviations, whereby the observed colour is much bluer, in some epochs, compared to model expectations (see Fig. 8 caption). The spectral index in these faint data points, instead of decreasing to a value of α ∼ −1 at g > 19 mag, diverges away from the model, to values of α = 0-0.3. A spectral index of +1/3 is expected for the overlapping radii of a viscously heated disk (eg. Frank et al. 2002). It may be that reprocessing on the disk becomes less important at these lower luminosities, revealing the viscously heated disk as the outburst fades.
If the viscous disk, with α = +1/3, is responsible for the deviations from the model, one would expect this to affect the shorter wavelengths more than the longer wavelengths, since this component rises at shorter wavelengths. This seems to be the case, with a colour deviation of ∼ 0.7-1.0 mag in the U -g and U -i CMDs, and ∼ 0.5 mag in the g -i CMD. The companion star could start to contribute to the optical emission at low fluxes, but we consider this to be unlikely to cause the observed deviations because (a) the star would have to be rising towards the blue, requiring it to be a more massive companion than is likely in this LMXB, and (b) the fluxes during the decay are still a couple of magnitudes above the quiescent level (see below), so the star is unlikely to dominate the emission. Optical emission from the viscous disk also has a shallower relation with the Xray flux, compared to reprocessing, and this is hinted at in Fig. 7 whereby the correlation slopes seem to appear shallower at lower luminosities compared to higher luminosities.
In the g -i CMD (see Fig. 8c), some of the brightest epochs show data that deviate from the blackbody model in the opposite sense; some data points are redder than the blackbody model by up to 0.2 mag in colour. This is less prominent in the other CMD filter combinations. Since g -i is the combination with the longest wavelengths, this is likely due to an additional component that is redder than the disk component, and which only makes a contribution at high luminosities. It is also variable; some data points are close to the blackbody model and so this redder component seems to vary in time. This is therefore probably the jet making a weak contribution to the i -band, since we know that the jet makes a stronger contribution at longer wavelengths in the infrared (Section 3.3) and it is variable (Section 3.2).

Long-term monitoring and quiescent magnitude
The only report of any optical quiescent magnitudes of GRS 1716−249 in the literature is a single weak constraint of B∼ 21.0-21.5 mag (della Valle et al. 1993). The source was not detected with Gaia in quiescence, and it only appears in EDR3 after the new outburst data were included (see Section 2.1.2). This provides a 20.7 mag limit in the G-band (Brown et al. 2016).
We have been monitoring GRS 1716−249 in quiescence with LCO (mostly using the 2-m Faulkes Telescope South) for the last 15 years, since 2006 February 3 (MJD 53769; see section 2.1.1 for details). The monitoring continues past the data we report on here through 2022 February. During quiescence, all the measurements obtained with XB-NEWS are forced-photometry points centred at the position of GRS 1716−249. On visual inspection of the quiescent data, we find that in most of the images the target is not visible at its expected position, and the quiescent magnitudes from XB-NEWS could be contaminated by emission from a brighter source very close to (∼ 2 arcseconds away from) the transient within the aperture, and a faint star 1.6 arcseconds away from the X-ray binary (see Fig. 1), making them unreliable.
To obtain a reliable quiescent optical magnitude, we select all the LCO images with good seeing (< 1.6 arcseconds) and perform aperture photomertry at the source position using an aperture size of ∼1 arcsecond in order to exclude the flux from the nearby stars and obtain the quiescent magnitude. The finding chart in quiescence obtained from the image with the best seeing (∼0.82 arcseconds) is shown in the right panel of Fig. 1. We find uncontaminated detections of the source at 13 epochs during quiescence spanning a range of 10 years (2011 May -2021 May, see Fig. 9). There is a slight variation during quiescence, with the i -band magnitudes ranging from 21.04±0.17 (MJD 59137.4) to 21.88±0.28 (MJD 58344.5). By combining all the 7 detections with seeing < 1.1 arcseconds, we find the quiescent i -band magnitude = 21.39±0.15 mag. The position of the source in quiescence is consistent with that measured from outburst data. The quiescent magnitude does not appear to change (within errors of the variations) before versus after the 2016-2017 outburst. We also rule out there being any long mini-outbursts after the 2017 outburst.
In the near-IR wavelengths, the source is not detected with 2MASS during quiescence (Rout et al. 2021), inferring an upper limit of 15.8 mag for the J band, 15.1 mag for the H band and 14.3 for the K band (Skrutskie 2006). From archival J-band images of the field taken during quiescence on 1999 July 5 and 7 (MJD 51364 and 51366) with the SOFI instrument at the New Technology Telescope (NTT; La Silla, Chile), we find a 3σ upper limit of the source of J > 18.8 mag, while the nearby southern star was found to have magnitudes of J = 15.38 ± 0.05 mag (see Section 2.2.2). There is a mention of probable near-IR quiescent magnitudes of GRS 1716−249, as Chaty et al. (2002) detected the source in J, H and K bands with the 2.2-m La Silla Telescope (ESO, Chile) on 1997 July 19 (MJD 50648) when it was expected to be in quiescence. They tabulate the quiescent magnitudes of the source as J = 19.4 ± 1.2; H = 19.2 ± 1.0; K = 18.3 ± 1.0. But they note that the source was not detected on 1998 July 6 (MJD 51000), and caution that observations with more powerful telescopes are needed to confirm the quiescent magnitudes.

DISCUSSION
Compared to the origins of X-ray or radio emission in a BHXB, the origin of the optical and near-IR emission is much less understood. Many physical processes could potentially contribute to the emission at these wavelengths, including X-ray reprocessing by the outer accretion disk (Cunningham 1976;Vrtilek et al. 1990), intrinsic thermal emission from a viscously heated outer accretion disk (Shakura & Sunyaev 1973;Frank et al. 2002), synchrotron emission originating from a steady compact jet during the hard state (e.g. Markoff, Falcke & Fender 2001;Jain et al. 2001;Corbel & Fender 2002;Buxton & Bailyn 2004;Russell et al. 2006;Kalemci et al. 2013;Saikia et al. 2019) and sometimes in transitional states (e.g. Fender, Belloni & Gallo 2004;van der Horst et al. 2013;Koljonen et al. 2015;Russell et al. 2020), a hot inner flow during the hard state (e.g. Veledina, Poutanen & Vurm 2013), and the companion star during quiescence (e.g. Casares & Jonker 2014). In this section, we explore the various emission processes contributing to the optical/UV fluxes of GRS 1716−249 using information from the methods mentioned previously, and discuss their implications on the system parameters, especially the distance to the source.

Optical/UV/IR emission mechanism
We study the optical/UV as well as broadband SEDs of GRS 1716−249 with quasi-simultaneous (within 24 hours) data, and find that they show a flat spectrum at optical/UV wavelengths (with a slight peak in the optical), with a positive slope in the near-IR regime, suggesting that the optical/UV emission mainly originates from a multi-temperature accretion disk. The optical/UV emission is near the peak of the blackbody from reprocessing. The fainter near-IR emission compared to optical is consistent with the Rayleigh-Jeans tail of the blackbody from the outer disk. The mid-IR emission on one date is comparatively brighter than what is expected from the disk alone. Such excess of emission in the IR regime is seen in many BHXBs (e.g., XTE J1550−564, Jain et al. 2001;4U 1543−47, Buxton & Bailyn 2004H1743−322, Chaty et al. 2015; XTE J1650−500, Curran et al. 2012;GX 339−4, Corbel & Fender 2002, generally associated with a compact jet. Along with being above the disk model, the mid-IR emission is also highly variable (see Table 1), and the radio to mid-IR spectrum is slightly inverted (with an index of α=0.13±0.03), which are typical signs of jet emission from BHXBs. The broadband spectral fitting performed by Rout et al. (2021) also shows that an irradiated accretion disk dominates the ultraviolet and optical emission. They report an IR excess compared to what is predicted by the irradiated disk model, and interpret it as due to the presence of a jet. Similarly, Bassi et al. (2020) fitted their broadband SED with the irradiated disk model diskir to describe the contribution of the accretion flow emission, which accounts for the irradiation of the outer disk and the reprocessing of the X-ray photons in the optical/UV band.
The slope of the optical/X-ray correlation also reveals the dominant emission mechanism of the accretion disk.
For an X-ray reprocessing accretion disk the slope of the correlation is theoretically expected to be ∼0.5 (van Paradijs & McClintock 1994). But we note that the theoretical value can be slightly different if there are extra contributions coming from additional emission components like irradiation from a disk wind, and can have a much larger range of slopes depending on which wavelength is used and whether the optical emission is coming from the Rayleigh-Jeans tail (RJ) or closer to the peak of the blackbody disk (Tetarenko et al. 2020;Shahbaz et al. 2015;Coriat et al. 2009). On the other hand, for a viscously-heated disk the slope of the correlation is expected to have a wavelength-dependent value ∼0.3 (Russell et al. 2006), and for an optically thick jet the expected slope is ∼0.5-0.7 (Corbel et al. 2003;Russell et al. 2006). For GRS 1716−249, the best-fit powerlaw correlations indicate the optical/X-ray slope to be 0.41±0.03 (see Section 3.4). This value is consistent with the X-ray irradiated accretion disk (van Paradijs & McClintock 1994), with additional contribution from the viscous disk, which could lower the value of the fitted slope from the theoretical value of ∼ 0.5 for irradiation (Russell et al. 2006). An X-ray irradiated accretion disk with optical emission coming from the peak of the blackbody, is also favoured by recent studies that find that the expected slope in the hard state can range from 0.13 (optical flux at RJ tail) to 0.33 (flux in the multicolour disk blackbody) for a viscously heated disc, and from 0.14 (RJ tail) to 0.67 (disk) for X-ray reprocessing with an isothermal disk (for a detailed calculation, see Tetarenko et al. 2020;Coriat et al. 2009). For cases like GRS 1716−249, where the outer disk temperature rises to ∼10000 K in outburst (see Fig. 8), the optical flux is found at the spectral transition between the RJ tail and the multicolour blackbody (Russell et al. 2006), and hence the optical/X-ray slope of 0.41±0.03 is consistent with the scenario of X-ray reprocessing. Similar values of power-law correlations have also been seen in other X-ray binaries like XTE J1817-330 (0.47± 0.03, Rykoff et al. 2007), GX 339-4 (0.44± 0.01, Coriat et al. 2009), GS 1354-64 (∼0.4-0.5, Koljonen et al. 2016). On the other hand, many sources like Swift J1357.2-0933 (Armas Padilla et al. 2013), Swift J1910.2-0546 , SAX J1808.4-3658 (Patruno et al. 2016) and Cen X-4 (Baglio et al. 2022) show a significantly shallower correlation (∼0.1-0.3). A slightly steeper correlation (∼0.56) is seen for V404 Cyg (Bernardini et al. 2016;Hynes et al. 2019;Oates et al. 2019), probably arising from contamination in optical fluxes from jet contribution. From our multi-wavelength correlation and spectral energy distribution analysis, we can rule out a significant optical emission component arising from a jet, in GRS 1716−249. This is also supported by our variability studies. Generally, sources with strong optical/IR variability on short (seconds-minute) timescales are known to have a strong jet contribution, and the variability is stronger at longer wavelengths where the disk makes a smaller contribution (Gandhi 2009;Gandhi et al. 2010;Baglio et al. 2018;Tetarenko et al. 2021). Disk variability is driven by changes in the mass accretion rate, which happen on the viscous timescale (days-weeks) for the viscously heated disk, and shorter (minute) timescales for reprocessing on the disk surface, if the X-rays have strong variability (with the reprocessing being smeared). The lack of strong variability in our optical data on short (minute) timescales, along with the presence of correlated variability on longer (days) timescales, suggests that the disk is producing the optical emission, and the contribution of the synchrotron jet emission at optical wavelengths is low in GRS 1716−249. The emission at near-IR wavelengths is dominated by the accretion disk, with a weak and variable jet component contributing towards the K-band in a few epochs. At mid-IR wavelengths, we find evidence for a highly variable jet component as suggested by the variable emission and the mid-IR to radio spectral index.
In addition to this, we investigated the color evolution of the source during its outburst. Our CMD analysis shows that the observed optical data mostly agrees with the single temperature blackbody model (at least at higher luminosities), with a scatter of ±0.1 mag in colour. This agreement supports the finding that the optical emission is originating mainly from a disk with varying temperature. The disk temperature varied between ∼7,000 K to ∼12,000 K, which is optimal for ioniz- ing hydrogen in the disk. At the brighter epochs in the CMD filter combination with the longest wavelengths (g -i ), we found that the data are slightly redder and brighter than what is expected from the disk model, with possibly some contribution coming from a jet. At the fainter epochs, we found significant deviations of the data from the reprocessing model, where the observed colour was much bluer, shifting the spectral index from α = 0 to +0.3 (which is expected in a viscously heated disk, Frank et al. 2002).
It is also worth considering the possibility of optical or IR emission from the hot flow. In this scenario, synchrotron emission from overlapping components of the hot flow contribute to the optical emission (Veledina, Poutanen & Vurm 2013). We have found that the optical spectrum is well described by a multi-temperaure disk, with the RJ tail in the IR (Section 3.3). The irradiation peak is detected, and there is low short term variability (Section 3.2). These characteristics, along with the behaviour in the optical/X-ray correlation and the CMDs, strongly favor a disk origin. In addition, Cuneo et al. (2020) detected variable, double peaked emission lines from hydrogen and helium, in the optical spectrum. These lines originate in the rotating accretion disk, and P-Cygni profiles were also detected from a disc wind. The hot flow model predicts strong short timescale optical variability (stronger in the optical compared to the IR because variations are amplified closer to the black hole, and the IR synchrotron emission in the hot flow orginates at larger radii) and a flat optical spectrum with a downturn at longer wavelengths (Veledina, Poutanen & Vurm 2013). We find stronger, high amplitude variability in the IR compared to optical, with the mid-IR flux density being higher than the optical on one date. The spectrum, emission lines, and variability properties are therefore not consistent with expectations from the hot flow. The hot flow model predicts a lower flux in the IR because the synchrotron-emitting region is physically limited by the inner edge of the accretion disc. So, while the optical is dominated by the disc, the mid-IR must be dominated by the jet, not the hot flow.

Constraints on the system parameters
We also conduct a comparative study of the quasisimultaneous optical/X-ray emission of GRS 1716−249 against a large sample of black hole and neutron star LMXBs in Fig. 10, with data taken from Russell et al. (2006Russell et al. ( , 2007. Both these classes of LMXBs are known to show different correlations, with the neutron star LMXBs being around 20 times optically fainter than black hole LMXBs (Russell et al. 2006), for reasons discussed in Bernardini et al. (2016). Assuming the previously estimated distance of 2.4±0.4 kpc (della Valle et al. 1994) is correct, GRS 1716−249 is found to be much more optically faint (or X-ray bright) compared to other BHXB samples (see the left panel in Fig. 10); in fact, at this distance GRS 1716−249 agrees more with the neutron star track in the global optical/X-ray correlation plot.

BH nature of GRS 1716−249
A BH nature of the source was first inferred by Masetti et al. (1996), who derived a lower limit for the compact object mass of >4.9 M from the superhump period of 14.7 hrs. Superhumps generally appear in disks with viscous-dominated emission, where the luminosity variations are caused by viscous dissipation associated with tidal deformation of the disk when it reaches the 3:1 resonance radius. We have shown in Section 4.1 with several lines of reasoning that the disk emission in GRS 1716−249 during outburst is dominated by X-ray irradiation. Such systems can have orbital modulations due to irradiation, rather than (or in addition to) superhumps (see the discussion in Haswell et al. 2001), especially when part of the optical variability comes from the irradiated face of the donor star. The superhump variability can be dominant at high orbital inclinations when the donor star to BH mass ratio is low and the donor star is shielded from irradiation, but optical modulation can be expected when the ratio is higher (see the discussion in Torres et al. 2021). As it is not clearly known whether the optical variability reported by Masetti et al. (1996) was a superhump or an irradiation effect, we cannot use it as a reliable constraint to the BH mass. Masetti et al. (1996) also noted that a massive primary is expected from the very long decay time of the X-ray light curve. Later, Tao et al. (2019) studied three quasi-simultaneous NuSTAR and Swift datasets of the system in its hard intermediate state, and assuming a distance of 2.4 kpc, constrained the upper limit for the compact object mass to be <8.0 M , at a 90% confidence level. Chatterjee et al. (2021) also used X-ray spectral analysis of the source during outburst to suggest a BH nature of the compact object. They fitted the X-ray spectra of GRS 1716−249 with the physical two-component advective flow (TCAF) model, keeping the mass of the primary as a free parameter, and constrained the mass of the compact object in the range of 4.5-5.9 M , but the values obtained are highly modeldependent and very unlikely to be a realistic range. The lack of Type I bursts during the outburst despite the presence of hydrogen (as suggested by the Hα lines, Cuneo et al. 2020) also provide strong evidence against a NS accretor. Moreover, Tao et al. (2019) show that good quality NuSTAR X-ray spectra of the source in the intermediate states can be fitted by BH models. In addition to the previous evidence, the X-ray timing properties of GRS 1716−249 also suggest that the compact object of the system is a BH. Chatterjee et al. (2021) report different power density spectra (PDS) of GRS 1716−249 in their Fig. 4, all of which show a strong decline from ∼3 Hz to 10 Hz. This behaviour is more typical of BH systems since the PDS of BHXBs show a strong decline at frequencies above 10-50 Hz (Sunyaev & Revnivtsev 2000). NS systems, on the contrary, can show variability up to 500-1000 Hz. The lack of X-ray pulsations and kilo-herz QPOs in the PDS (typical signatures of neutron star systems), the presence of type-C and type-B QPOs in the PDS of GRS 1716−249 (Chatterjee et al. 2021), and the strong decline of the power spectra below 10 Hz all reinforce the identification of the compact object as a BH.
If GRS 1716−249 is indeed a black hole, then the discrepancy shown by GRS 1716−249 with respect to other BHXBs in the global optical/X-ray correlation space could have two possible explanations. Either the source is intrinsically much more optically faint than what has been observed in other BHXBs at a given X-ray luminosity, or it is located much further away than was previously thought.

Distance to GRS 1716−249
The original distance estimate of 2.4±0.4 kpc is based on a comparison of the source to other X-ray binaries with data from a few decades before (della Valle et al. 1994). They argued that the lower limit on the distance is expected to be ∼2 kpc from the equivalent width of the NaD absorption lines. To constrain the upper limit, the peak optical brightness was compared to other BHXB outbursts known at the time (della Valle et al. 1994). Since then, a distance of 2.4±0.4 kpc has been used by various studies concerning GRS 1716−249. Later Hynes (2005) notes that one should be cautious about using such method to constrain the upper limit on the distance. It has now become clear that BHXB outbursts can peak at different luminosities, from close to the Eddington limit, down to ∼ 10 36 erg s −1 or less (e.g. the Very Faint X-ray Binaries, or mini-outbursts; Heinke et al. 2015;Zhang et al. 2019). Moreover, the historic peak optical brightness of LMXBs during outburst used by della Valle et al. (1994) was based on a compilation that neither corrected for orbital period, nor performed sorting of neutron stars versus black holes (van Paradijs 1981). In light of all these arguments we do not consider the formerly estimated upper limit of 2.8 kpc as a reliable constraint. Masetti et al. (1996) had discovered optical modulations in the source with a prominent period of ∼ 14.7 hrs, and found that the secondary star in the system should be substantially brighter than claimed by della Figure 11.
Expected peculiar velocity (υpec) of GRS 1716−249 for a range of possible distances over 0.5-30 kpc. The shaded region bound by the two dashed lines represent the 1σ scatter. The kinematics of the system favour distances of 8 kpc, as for BH systems the natal kick probability distribution ranges only up to ∼150 km/s (Mandel & Muller 2020). Valle et al. (1994). To explain this discrepancy, they suggested that either the distance of 2.4±0.4 kpc has been underestimated, or the secondary is a slightly evolved late-type star.
In addition to the previous arguments, we also find that a higher value of distance is expected from the state transition luminosity distribution of the source (e.g. Maccarone 2003;Kalemci et al. 2013;Vahdat Motlagh et al. 2019). It has been observed that BHXBs transit from the soft state to the hard state at luminosities between 0.3-3% percent of the Eddington luminosity (Kalemci et al. 2013), with a mean value of 1.9±0.2% (Maccarone 2003). The state transition luminosity has been used to estimate the distances to many BHXB sources (e.g. Homan et al. 2006;Miller-Jones et al. 2012). Although GRS 1716−249 did not go to a soft state, we use the luminosity during transition from the final hard/intermediate state to the hard state (MJD 57978, Swift/XRT flux 9.25 × 10 −10 erg s −1 cm −2 at 2-10 keV) to estimate the distance. Assuming a BH of mass 7M , 1.9±0.2% Eddington luminosity and a bolometric correction factor of 2 relative to the Swift/XRT band, we obtain a probable distance of 8.7±0.5 kpc for the source. For a more conservative range of 0.3-3% percent Eddington luminosity (Kalemci et al. 2013), the distance range is found to be 3.46-10.94 kpc.
Moreover, for a distance of 2.4 kpc, the inner disk radius depending on the inclination angle is r in ∼ 15 km (see Fig. 6 of Bassi et al. 2019), which is very unusual for a BH disk spectrum; while a more plausible value of r in > 50 km is obtained for distances d > 8 kpc. An underestimated distance could also explain the discrepancy we see for this source with respect to other BHXBs in the optical/X-ray correlation plots (see Fig. 10). From our global correlation comparison, we find that for a distance of 4 kpc and less the data are more consistent with being a neutron star, and for distances more than 4 kpc the data are more consistent with a BH.
We also place a conservative upper limit on the distance as 17 kpc from the global optical/X-ray corelation plot (see Fig. 10), as for a greater distance, the source would be the most X-ray luminous BHXB, probably exceeding the Eddington limit (depending on the black hole mass). The proper motion estimate of the source is ∼ 4.65 ± 1.12 mas/year, and the potential kick velocity (after removing Galactic rotation) is ∼70-100 km/s for a distance d =2.4 kpc (Atri et al. 2019). We performed a simulation using all the standard assumptions of Gandhi et al. (2019) and the measured proper motions assuming a radial velocity of -10 km/s (for further details, see Atri et al. 2019), and found that at any distance higher than 6 kpc, the space velocity of the source starts to exceed 100 km/s. At the Galactic centre distance, the source peculiar velocity increases to ∼150 km/s (see Fig. 11). On the other side of the Galaxy, however, median peculiar velocities are predicted to be between ∼190-330 km/s. Such high velocities are not expected in BH systems, where the natal kick probability distribution ranges up to ∼150 km/s with a root-mean-square kick of ∼60 km/s (Mandel & Muller 2020). So the kinematics of the system favour distances of 8 kpc, suggesting that it is significantly closer than 17 kpc. We also note that the source had a failed-transition outburst and did not show a transition to the soft state (Bassi et al. 2019). Generally, the failed-transition outbursts reach lower peak X-ray luminosities than full outbursts (Tetarenko et al. 2016;Alabarta et al. 2021). For example, in the case of one of the best studied BHXB GX 339−4, the luminosity at which the hard-to-soft state transition occurs during a full outburst is ∼ 0.11L Edd , and the luminosity during failed-transition outbursts are always equal to or lower than this value (Tetarenko et al. 2016). Assuming a similar behaviour of ∼10% L Edd during the peak flux in the case of GRS 1716−249, a conservative mass of 7M results in a distance of 5.0 kpc, and adds additional support to a closer distance. From all the arguments stated above, we place a conservative upper limit of 17 kpc for the system, although our lines of evidence suggest a much lower value (∼8 kpc).
From the global correlation plot and the list of reasoning mentioned, we constrain the distance of  Figure 12. Radio/X-ray luminosity correlation plot of GRS 1716−249 with quasi-simultaneous data from Bassi et al. (2019), and a large sample of black hole and neutron star LMXBs with data made available by Bahramian et al. (2018). The X-ray luminosities are in the 1-10 keV energy range, while the radio luminosities are taken in ∼5 and 9 GHz. The green crosses depict the position of the source assuming the distance of 2.4 kpc (literature value), and the red stars represent the position when a greater distance (in this case, 8 kpc) is assumed.
GRS 1716−249 to likely be in the range of 4-17 kpc (see the right panel in Fig. 10), with a most likely range of ∼ 4-8 kpc. This improved distance estimate will have implications for models of GRS 1716−249 that depend on its distance (e.g. Bassi et al. 2020;Zhang et al. 2021;Chatterjee et al. 2021); and affect the parameters of GRS 1716−249, like the inferred masses, spins and inclination angles, which depends critically on the assumption of d =2.4 kpc (Tao et al. 2019). For example, the same spectral modelling assuming a distance of d =8 kpc results in M BH = 24.8 +1.7 −10.2 M (compared to M BH = 7.6 +0.8 −2.7 M obtained assuming d =2.4 kpc), which shows that the parameters of the system are clearly model dependent and sensitive to the distance (Lian Tao, private communication). So distances beyond 8 kpc would make it the most massive stellarmass BH known in our Galaxy (Cyg X-1 currently holds the record with 21.2±2.2 M , Miller-Jones et al. 2021), which is unlikely. This implies that a distance above 8 kpc is highly implausible and again argues for an upper limit that is substantially lower than 17 kpc. A greater distance than the literature value also changes the position of the source in the radio/X-ray luminosity correlation plot (Bassi et al. 2019). Generally, BHXBs follow two different tracks in the radio/X-ray luminosity correlation plot (where L R ∝ L β X ) -the 'standard' track with a power law index β ∼0.5-0.7 (e.g. Corbel et al. 2003Corbel et al. , 2013Gallo, Degenaar & van den Eijnden 2018), and the much steeper 'outlier' track with β ≥1 (e.g. Corbel et al. 2004;Coriat et al. 2011;Gallo, Miller & Fender 2012), although the existence of two separate tracks has been questioned statistically (e.g. Gallo et al. 2014;Gallo, Degenaar & van den Eijnden 2018). Newer studies have also shown that the two tracks are not welldefined and there is evidence for standard track sources to get steeper at high X-ray luminosity (Koljonen & Russell 2019), and outlier track sources to get shallower and ultimately rejoin the standard track at low X-ray luminosities (see e.g. Coriat et al. 2011;Carotenuto et al. 2021). Although the underlying physics behind the two tracks is not completely clear, it has been suggested that the dichotomy could originate either from the structure of the inner accretion flow or from different physical properties of the jets resulting in different levels of radio emission (Coriat et al. 2011). NSXBs, especially in the hard state, also showed a similar correlation with fainter radio emission compared to BHs and a steeper slope (with β ∼1.4, e.g. Fender & Hendry 2000;Fender & Kuulkers 2001), although there is strong evidence of different classes of NSXBs showing different behaviour in the radio/X-ray correlation plane (see e.g. Tudor et al. 2017;van den Eijnden et al. 2021). We find that a greater distance shifts GRS 1716−249 from the lower luminosity part of the outlier track (luminosities where mostly NSs are observed and the BHs seem to shift to the standard track) to the higher luminosities where the majority of the BHXB sample following the outlier track lies (see Fig. 12). A greater distance also implies that the source was potentially formed in the bulge, and hence its proper motion is not necessarily representative of its natal kick, since the bulge itself has a large velocity dispersion and scale height (Atri et al. 2019).

CONCLUSION
The 2016-2017 outburst of the BHXB GRS 1716−249 (or GRO J1719−24) is well-studied at X-ray and radio wavelengths. In this work, we investigate the optical, near-IR, mid-IR and UV wavelength monitoring data of GRS 1716−249 in outburst using LCO, REM, VLT (VISIR) and Swift's UVOT, and compare them with the multi-wavelength archival data from Gaia, Mount Abu 1.2 meter telescope, Swift XRT, NuSTAR, MAXI, ATCA, VLA and LBA. We also report the long-term (∼ 10 years) optical light curve of the source using LCO and find that the quiescent i -band magnitude is 21.39±0.15 mag.
We find that the optical and UV emission of the source in outburst is mainly originating from a multitemperature accretion disk, with X-ray reprocessing dominating at high luminosities, and with some con-tribution at the fainter end from the viscously heated disk. Although the near-IR emission is dominated by the emission from the accretion disk, it has a weak contribution from the variable jet in a few epochs in the K-band. The mid-IR and radio emission of the source are dominated by the synchrotron emission from a compact jet. In the hard state, the optical/UV emission of the source is correlated with both the soft and hard X-ray emission. The power-law coefficient of the correlation is consistent with the optical emission coming from an X-ray irradiated accretion disk with possibly some additional contribution from the viscous disk, as a hint of a shallower coefficient at low luminosities. This is also supported by the spectral energy distributions, variability studies and color-magnitude diagrams of the source during the outburst.
Finally, we discuss how the previous estimates of system parameters of the source (especially its mass and distance) are based on various assumptions, and cannot be completely trusted. From the global optical/X-ray correlation study in comparison with other black hole and neutron star X-ray binaries, and several other lines of reasoning, we show that GRS 1716−249 is much further away than what has previously been assumed, with a probable distance within the range 4-17 kpc, and a most likely range of ∼ 4-8 kpc.

ACKNOWLEDGMENTS
The authors thank the anonymous referee for useful comments and suggestions. The authors also thank Lian Tao