© ESO 2021

1 Introduction

High-contrast imaging with extreme adaptive optics (ExAO) facilities at JHK-band wavelengths has yielded many sub-stellar companions to date, including some directly imaged planets (e.g. Marois et al. 2010; Lagrange et al. 2010; Macintosh et al. 2015; Keppler et al. 2018). Since giant exoplanets are fairly hot (~10³ K) at young ages (~10⁷ yr) after their formation (Spiegel & Burrows 2012; Marleau et al. 2019), they emit a large fraction of their bolometric flux at near-infrared (NIR) wavelengths. This fact, along with the relatively low thermal background in the JHK range and the fairly high adaptive optics (AO) correction quality that can be reached there, is the reason most exoplanets searches – and most detections – occur at JHK wavelengths. However, most exoplanets imaged thus far are relatively young, with an age of ~1–100 Myr. Older gas giant planets, at approximately gigayear ages, will have typical effective temperatures of ~100–300 K. At such temperatures, the flux at JHK wavelengths becomes so low that it is out of reach for present-day facilities. For a range of atmospheric properties, a significant flux bump remains at 4 μm (e.g. Allard et al. 2001; Burrows et al. 2006; Fortney et al. 2008), which is included in the red end of the L′ filter. Hence, L′ is the shortest feasible wavelength band where old and cold planets can be studied. At longer wavelengths, the planet-to-star contrast is often even more favourable, although the thermal background also increases rapidly with wavelength longwards of L′, posing a considerable observational challenge of its own. As a result, a number of efforts have been made to develop high-contrast imaging for detecting and characterising planets in L′ from the ground (e.g. Kasper et al. 2007; Janson et al. 2010; Quanz et al. 2010; Absil et al. 2013) and in space (e.g. Janson et al. 2015; Durkan et al. 2016; Baron et al. 2018).

An important goal in exoplanet research is to detect planets both in imaging and with radial velocity (RV) and/or astrometry. This allows us to simultaneously determine a range of properties such as the mass, orbit, luminosity, and spectral distribution of the planet, providing a much larger information space than is available with either technique in isolation. For brown dwarfs, a combination of imaging and RV (e.g. Crepp et al. 2012, 2014; Peretti et al. 2019) or imaging and astrometry (e.g. Calissendorff & Janson 2018) or all three (e.g. Brandt et al. 2019; Grandjean et al. 2019; Maire et al. 2020b; Currie et al. 2020) has been achieved in a number of cases, encouraging increasing efforts in recent years to use similar approaches to detect and characterise exoplanets (e.g. Mawet et al. 2019). In two cases so far, directly imaged planets have also been observed astrometrically. The planet β Pic b (Lagrange et al. 2010) causes an acceleration of its host star between the Hipparcos (Perryman et al. 1997) and Gaia (Gaia Collaboration 2018) astrometric epochs (Snellen & Brown 2018; Nielsen et al. 2020; Brandt et al. 2021). In combination with direct monitoring of the planet motion, this yields a dynamical mass estimate of 11 ± 2 M_J. Additionally, Nowak et al. (2020) recently combined direct observation of a second planet in the system, β Pic c (Lagrange et al. 2019), with RV and astrometric data, constraining its inclination and luminosity to estimate its mass at 8.2 ± 0.8 M_J. β Pic is a young (~24 Myr; see e.g. Bell et al. 2015) system and very unusual in the sense that it is also very nearby at 19.44 pc (Gaia Collaboration 2018) – most similarly young planetary systems are much more distant and thus harder to detect astrometrically. This is further emphasised by the fact that they can only be imaged at very large separations, requiring excessively large astrometric baselines for dynamical detection. The most promising overlaps between the dynamical and imaging characterisation of planets is in very nearby systems; however, such systems are generally relatively old and thus their planets can be expected to be relatively cold, further motivating the need for direct imaging developments in the L′ band and longwards (Heinze et al. 2010).

Located only 3.639 ± 0.003 pc away (Gaia Collaboration 2018), ϵ Ind A is one of the most nearby Sun-like (K5-type, 0.76 M_sun) stars (e.g. Demory et al. 2009). A common proper motion low-mass companion, ϵ Ind B, was reported at a wide separation (~1460 au) in 2003 (Scholz et al. 2003), which shortly afterwards was discovered to itself be binary, forming the brown dwarf pair ϵ Ind Ba and Bb (McCaughrean et al. 2004). Age estimates in the literature for ϵ Ind A from chromospheric activity and rotation have ranged from ~1 to ~5 Gyr (e.g. Lachaume et al. 1999; Barnes 2007; Feng et al. 2019). However, the lower age estimate of 1 Gyr by Barnes (2007) is based on an older, less accurate value of the rotation period for the star, compared to the more recent and reliable estimate of ${35.732}_{-0.003}^{+0.006}$ days by Feng et al. (2019), which points to an older age of ~4.5 Gyr using the same method of gyrochronology (Delorme et al. 2011) as used by Barnes (2007). This indicates that the age of ϵ Ind A system is expected to be older. Additionally, ϵ Ind A has been observed to exhibit a long-term RV trend, which was stronger than the secular acceleration caused by the rapid motion of the star on the sky (Endl et al. 2002; Zechmeister et al. 2013). The trend was also too strong to be explained by any dynamical impact of ϵ Ind Ba/Bb, indicating that an additional object – most likely a giant planet – must be exerting a small but significant gravitational pull on the star. Hence, dedicated imaging surveys were performed (Geißler et al. 2007; Janson et al. 2009) in order to try to detect the companion but yielded no detections. This further underlined that the companion must be a planet since the detection limits excluded brown dwarfs and more massive objects as potential companions. Meanwhile, RV monitoring continued for the object and was showing increasingly clear curvature, underlining the reliability of the RV solution and greatly enhancing the predictability of orbit and mass for the planet. On this basis, deeper imaging campaigns were planned and executed, which are the topic of thisarticle. Recently, a sufficient coverage in both RV and astrometry was reached such that the planet ϵ Ind Ab could be confirmed and all of its orbital elements constrained (Feng et al. 2019). The current best-fit mass and semi-major axis of ϵ Ind Ab are 3.25${}_{-0.65}^{+0.39}$ M_J (at an inclination of ${64.25}^{°}{}_{-6.09}^{+13.8}$) and 11.55${}_{-0.86}^{+0.98}$ au.

Here, we present deep AO imaging observations with the 3.8 μm NaCo (nasmyth adaptive optics system (NAOS) + near-infrared imager and spectrograph (CONICA)) L′ and 10–12.5 μm NEAR (new Earthsin the α Centauri region) wavelength filters. NaCo (Lenzen et al. 2003; Rousset et al. 2003), which was a workhorse NIR instrument at the Very Large Telescope (VLT) for nearly two decades, was recently decommissioned, and a new instrument with enhanced 1–5 μm imaging capabilities, ERIS (Riccardi et al. 2018; Dolci et al. 2020), will see its first light at VLT as early as this year. NEAR (Käufl et al. 2018; Pathak et al. 2021) was an upgrade of VISIR (Lagage et al. 2004), with AO assistance and a coronagraph added to the mid-infrared camera. NEAR was developed primarily for the purpose of searching for planets in the α Cen system (Kasper et al. 2019), but was briefly offered for broader scientific applications in a science verification programme executed during the second half of 2019 before also being decommissioned. Between them, the two instruments offer unprecedented contrast and sensitivity in a wavelength range in which cool-to-temperate planets radiate a significant portion of their thermal flux. We describe the observational setups in Sect. 2 and the data reduction schemes in Sect. 3. The results are discussed in Sect. 4, and the conclusions are summarised in Sect. 5.

2 Observations

2.1 NaCo observations

Observations with NaCo of ϵ Ind A were acquired in October and November of 2018, using the L′ filter in pupil tracking mode. Originally allocated in service mode, the observations were transferred to designated visitor mode when the visible wavefront sensor (WFS) broke down and a non-standard implementation of the infrared WFS was required due to the high brightness of ϵ Ind A. The first ~quarter-night run was executed on October 5, but conditions were poor and the data were not usable. This paper will focus on the subsequent four ~quarter-night runs performed on October 12 (MJD = 58 403), October 26 (MJD = 58 417), November 3 (MJD = 58 425), and November 4 (MJD = 58 426), all of which had a seeing of ~1′′ or less with clear sky conditions.

During the science sequence, 1 minute exposures were performed through single readout blocks of 300 frames (NDIT) already coadded at the detector level, each with an integration time (DIT) of 0.2 s. Full-frame readouts were used so that the star could be placed in the centre of a ‘good’ quadrant of the detector, avoiding the lower left quadrant which was affected by a large number of unusable pixel stripes. At the time of the observations, the ephemerides of ϵ Ind Ab were only loosely constrained, theoretically allowing for separations of up to ~5′′. Placing the star near the centre of the upper left quadrant of the ~27-by-27′′ detector left enough space both relative to the bad quadrant and the detector edges to accommodate detection at such large separations. The expected separation of ϵ Ind Ab during the NaCo observation epoch (i.e. as on October 26, 2018) was later determined to be 1.37′′ ± 0.16′′, following the ephemeridescalculated by Feng et al. (2019). The field rotation that was obtained for the data taken on October 12, October 26, November 3, and November 4 were 29.5°, 34.5°, 58.5°, and 58°, respectively.In total across the four usable observing runs, 456 science frames were acquired, of which 8 were visually classified as being of poor quality and discarded. This left 448 co-added frames for the analysis, corresponding to a total useful integration time of 3.73 h. No dithering was performed during the observations, but individual sky frames were interspersed among the science observations to monitor the thermal background. The observations were obtained in standard (saturated) imaging mode with pupil tracking. The resulting science data were saturated, typically out to ~ 6 pixels (0.16′′) from the central star. Non-saturated frames of the primary star were acquired with an integration time of 0.1s and 100 coadds, with the use of the ‘long’ neutral density (ND) filter of NaCo, with a transmission in the L′ -band of 1.8%¹.

2.2 NEAR observations

ϵ Ind A was observed using NEAR in visitor mode in September 2019, utilising the enhanced sensitivity and point spread function (PSF) contrast obtained via its combination of AO with an AGPM coronagraph (Maire et al. 2020a). Observations were carried out during three nights on September 14 (MJD = 58 740), September 15 (MJD = 58 741), and September 17 (MJD = 58 743) using the NEAR (10–12.5 μm) filter in pupil-stabilised tracking mode in good weather conditions. Science exposures were performed with high-frequency deformable secondary mirror (DSM) chopping at a frequency of 8.33 Hz and a chopping amplitude of 4.5′′, in combination with nodding. We used an integration time of 60ms – resulting in two chop images per nod – which we averaged. The star was positioned behind the coronagraph near the centre of the detector with an effective field of view of ~ 20′′ × 20′′ in the small field (SF) mode. The expected separation of ϵ Ind Ab as on September 15, 2019, was 1.02′′ ± 0.19′′ (3.71 au) at a position of ⟨RA⟩ = 0.903′′ ± 0.204′′, ⟨Dec⟩ = −0.302′′ ± 0.352′′ relative to the star.

The science frames for the target were obtained from the difference image between the two nod positions, giving an on-target integration time of 48s for each nod beam position. A total of 246 such nod difference frames were obtained over the three nights, out of which 23 frames were of poor quality due to unfavourable seeing, AO or tracking malfunction and were discarded. A total of 223 science frames were thus used for further analysis, corresponding to a total integration time of 2.23 h. An image of ϵ Ind A not obscured by the coronagraph was obtained for analysis from the nod-difference frames by summing the off-centre chop beams at the two nod positions and averaging over all such frames. The integration time for this unobscured image of the central star obtained via the above method is the same as the on-target integration time of each nod difference frame in the coronagraphic images. Since the star does not have a very high photon flux rate in the mid-infrared relative to the background, no specific measures were taken during observation to prevent saturation.

3 Data reduction

3.1 NaCo reduction

NaCo was installed at the VLT in 2001, so it had been operational for approximately 17 yr when the observations were acquired in late 2018, and was decommissioned less than a year later in 2019. By the time of the observations, NaCo was experiencing several technical difficulties, not least with respect to the detector, which had a large number of bad pixels and a rather unstable flat field with variable stripes and similar patterns. This required particular considerations in the data reduction procedure, expanding on previous NaCo L′ processing (e.g. Janson et al. 2008). As we will see, with dedicated processing it was possible to reach a higher image depth than any other observation of the system, including previous NaCo imaging (Janson et al. 2009), although not by the samefactor that would have been expected if the detector status had been as good as in those previous epochs.

Sky flats² are essentially the only way to get reliable flat field in L′ with NaCo, and are acquired on a monthly basis at the telescope. Flat field frames acquired during the last years of NaCo operations exhibited variable features such as horizontal stripes occurring at some epochs but not others; quasi-vertical stripes that changed slightly but significantly in morphology between epochs; and diffuse circular or elongated features across varying parts of the detector field. For this reason, we took special care in identifying a well-matching flat to the observational data, visually examining every flat field with relevant settings that had been acquired throughout 2019. A flat field from November (which is close to the epoch of the observations) was found to be the best option, although it had to be filtered from an additive horizontal stripe pattern that occurred in individual frames. This was done row-by-row by calculating the median flux of 50 pixels near the edges of the frame in that row and subtracting the result from each pixel in the row. Each science frame was dark subtracted and flat field corrected. Bad pixels were identified through their deviations from a 5-pixel median box filtered version of the flat field. Unsharp masking was applied to every science frame through a median box filtering with a box size of 20 pixels that was subtracted from the image. This eliminates the low spatial frequencies that dominate the distribution of the thermal background flux, while having a negligible effect on point sources that have a full width at half maximum (FWHM) of ~4 pixels. The exact same procedure was applied to the non-saturated frames, such that any small flux losses from the procedure would be accurately represented when calibrating the contrast in the final images from the unsaturated PSF.

Post-processing of the science frames was carried out using a combination of two different techniques. In the regime close to the star (<3′′), we used an automated version of the regime switching model (RSM) algorithm (Dahlqvist et al. 2021), which operates based on the temporal evolution of pixel intensity in the de-rotated cubes of residual frames resulting from an optimum combination of different PSF subtraction techniques based on angular differential imaging (ADI; Marois et al. 2006). In a recent data challenge designed to evaluate the relative performance of different high-contrast algorithms, RSM scored very highly both in terms of high true detection rates and in terms of low false detection rates (Cantalloube et al. 2020). Iterating over radial distance, the RSM algorithm uses a two-state Markov chain on the resulting time-series to generate probabilities, for each pixel in an annulus, to be in two regimes; one where the pixel intensity is described purely by residual noise from quasi-static speckle and another where the pixel intensity is described by both residual noise and a model of the planetary signal.

The RSM algorithm (Dahlqvist et al. 2020) relies on the off-axis PSF, or a forward-modelled PSF (Dahlqvist et al. 2021) to account for signal self-subtraction resulting from the reference PSF subtraction, to model the planetary signal. The RSM detection map is then produced by time-averaging the probabilities of being in the planetary regime. In this work we used an enhanced version of the RSM algorithm to generate our final contrasts. This version optimises the different PSF subtraction techniques used in the algorithm based on minimisation of contrast; optimises the algorithm itself based on maximization of probability ratio of injected fake companion to the background noise; and then searches for the best set of PSF-subtraction techniques and observation sequences to generate the optimum contrast. This best set is selected via a bottom-up greedy selection algorithm. The algorithm iteratively adds the PSF-subtraction technique and observation sequence that maximises the incremental increase in the probability ratio of injected fake companion to the background noise, until no more incremental increase can be obtained. Consequently, the final detection map shown in Fig. 1a was generated using the local low-rank plus sparse plus Gaussian decomposition (LLSG; Gomez Gonzalez et al. 2016) for the first observation sequence, the annular principle component analysis (Annular PCA; Gomez Gonzalez et al. 2017) for the second, the locally optimised combination of images (LOCI; Lafrenière et al. 2007)for the third, and a combination of all three PSF-subtraction techniques for the fourth observation sequence³ (respectively, October 12th, October 26th, Novermber 3rd, and November 4th). No companions were detected in the final map.

Since the noise statistics in the map are unknown, a standard 5σ threshold cannot be defined to generate contrast curves. The contrast curves were instead obtained by injecting fake companions into the observation sequences with different flux values, at a range of azimuths and radial separations. The detection limit was defined as the flux at which the true positive rate (TPR) corresponding to the detection of the first false positive in the entire frame (< 3′′) is 50% (which is similar to a standard 5σ criterion, see Dahlqvist et al. 2021 for more details). To add further robustness to our results, we also produced a contrast curve corresponding to a TPR of 90%. In Figs. 2a,b, we show examples of the 50 and 90% criteria by injecting 11 fake companions into the observing sequence at the level of the calculated contrast curve, all at a mutual azimuthal angle, along the same radial distances at which the contrasts were estimated. The figures show the detected signals at a threshold equal to the brightest speckle in the empty probability map (0.0025), at a TPR of 50 and 90%, respectively. Out of the injected companions, four were retrieved for the 50% case and nine for the 90% case. The slight deviations from exactly 50 and 90% retrievals were due to statistical fluctuations; the numbers quoted here correspond to a single azimuth, while fake companions were injected at ten different azimuths per angular distance to compute the contrast curves.

As we mentioned above, the RSM algorithm was applied out to a separation of 3′′ ; in the background-limited regime outside of this range, RSM offers no advantages relative to conventional ADI-based PSF subtraction techniques (Dahlqvist et al. 2021). Furthermore, owing to the aforementioned detector noise and instability of the flat field correction, we have found that conventional techniques such as ordinary LOCI perform significantly worse in this regime than expected, given the integration time and predicted background level.

To mitigate this effect, we performed an alternative LOCI-based subtraction procedure. In this alternative implementation, we ignored the small drift of the stellar PSF that occurred during the observations, and operated on the images as if the star had been perfectly fixed with respect to the detector. In other words, LOCI was applied to the non-shifted individual images, and shifting was applied only after the reduction had finished. This procedure yielded a worse performance than the normal procedure in the contrast-limited regime because the unaccounted relative shifts in thestellar PSF resulted in a worse average PSF matching. However, the large benefit of the non-shifted implementation was in the (nominally) background-limited regime. When the images were shifted prior to optimisation and subtraction, pixels were shifted from their original locations, such that any systematic error that might have occurred at the pixel level gets spread across different parts of the resampled images, forming a quasi-Gaussian noise term in each new pixel with (in this case) a larger variance than the shot noise of the thermal background. By operating on the non-shifted frames, the LOCI subtraction has the chance to subtract out systematic noise effects on the detector level, potentially yielding a better performance outside of the contrast-limited central regime. A 5σ contrast curve as a function of separation was derived based on the standard deviation of samples within an annulus of a given separation, where each sample was the aperture sum of pixels around the sampling point with a diameter equal to the FWHM, 120 mas. Self-subtraction was accounted for by using the unsaturated stellar image as a PSF representation and calculating the flux loss at every location in the image based on the actual optimisation coefficients used in the modified LOCI subtraction. For the final output, we used the modified LOCI implementation from 3′′ and outwards, which is outside of the region used for the optimisation between 2.7′′ and 3′′. Hence, there was no systematic subtraction of any hypothetical companions from the optimisation in the region of interest, and therefore the only recurring source of self-subtraction came from partial PSF overlaps from each pair of companion signatures in any pair of mutually subtracted images. Due to the large separation of >3′′, such overlaps are small, and as a result self-subtraction is low in the part of the special LOCI output image used in this analysis. The non-shifted procedure did indeed perform very well in the outer ranges of the image, with a substantially better contrast curve than the shifted procedure in that regime, as can be seen from Fig. 3, which shows the comparison of contrast curves obtained using different reduction techniques used on NaCo L′ data in this work. The combined contrast curve from both RSM and the modified LOCI procedure is shown in Fig. 4, with the sensitivity limits from RSM shown for both 50 and 90% TPR.

As an alternative check, we also reduced NaCo data by using a combination of the classical LOCI, the modified LOCI, and TRAP algorithms. In the classical version of LOCI we shifted the images to a common central pixel location before PSF subtraction, each shift determined through cross-correlation, and performed optimisation in an annulus between 40 and 60 pixels of separation from the central star. Self-subtraction was accounted for in the same way as described above for the modified LOCI reduction. As before, in the background limited regime, we performed the modified LOCI implementation where the individual images are not shifted to a common centre before reduction to account for the effects from systematic detector noise. The resulting 5σ contrast curve was also derived similar to the modified LOCI reduction. In the regime very close to the star, the LOCI procedures described above exhibit substantial self-subtraction, so we ran the TRAP algorithm (Samland et al. 2021) on the data, with data from the different nights treated as one contiguous dataset. TRAP being a temporal optimisation algorithm, instead of shifting the actual images, the stellar drift in the images was incorporated in the temporal forward model constructed by the algorithm, with the fraction f of principal components set to 0.3. The resulting final output of TRAP showed an improvement in contrast at the very smallest separations (within ~0.6′′), as expected.Hence, we derived contrast curves based on three separation ranges in which three different algorithms provide an optimal result. Consistent with the previous reduction, no companions were detected in the reduced image thus obtained either. However, the contrast curves derived in this way were up to an order of magnitude less sensitive than the contrastcurve from RSM algorithm at both 90 and 50% TPR (refer Fig. 3). Thus for the reminder of this analysis, we continued with the contrast curves derived from RSM and modified LOCI algorithms (Fig. 4). Conversion into absolute magnitudes was based on WISE (Wright et al. 2010) photometry, where the W1 filter was used as a proxy for the L′-band (m_W1 = 2.9).

Fig. 1 Final detection probability map from NaCo (panel a) and the final reduced image from NEAR observations (panel b) of ϵ Ind A. The predicted position of the companion at the time of observation, along with the 1σ uncertainty in position, as inferred from the best-fit ephemerides, based on RV and astrometry data (Feng et al. 2019), is shown asan ellipse. The NaCo image was reduced using a combination of RSM and modified LOCI subtraction applied at suitable separations while the NEAR image was reduced using the circular profile subtraction technique. No companions were detected in either image.

Fig. 2 Detection maps from NaCo reduction using RSM at 50 (panel a) and 90% TPR (panel b). Out of the 11 injected fake companions along the shown axis, four have been retrieved for the 50% case and nine for the 90% case. The detected signals correspond to a threshold of 0.0025, which is equal to the brightest speckle in the empty probability map.

Fig. 3 Contrast curves resulting from the different reduction algorithms applied to NaCo L′ observationsin this work. The y-axis shows the 5σ contrast limits obtained from TRAP, the normal LOCI, and modified LOCI subtraction techniques (shown as blue, green, and red dotted line), and the equivalent flux ratio or sensitivity limits obtained using RSM algorithm at 50 and 90% TPR as yellow and black dotted lines, respectively.

Fig. 4 Composite contrast curve from NaCo L′ observations, obtained from the combination of RSM and the modified LOCI subtraction techniques. The y-axis shows the 5σ contrast limits obtained from modified LOCI subtraction techniques (shown in red solid dots), beyond ~11 au and the equivalent flux ratio or sensitivity limits obtained using RSM algorithm corresponding to a TPR of 50% (blue solid dots) and 90% (green solid dots).

3.2 NEAR reduction

The science data from NEAR observations of ϵ Ind A, after binning, were available as chopping subtracted and nodding-corrected frames. Since chopping takes care of background as well as any bias or dark current in the data and nodding cleans up any residuals from change in the optical path through the telescope and instrument due to chopping, no additional measures were taken to correct for such effects. Following Petit dit de la Roche et al. (2020), no flat fielding corrections were applied to the data, as they do not improve the performance for NEAR data. The stellar PSF subtraction was thus performed on the 223 nod-difference frames stacked into three cubes, one for each observation night.

At the time of observation, the expected planet-star separation for ϵ Ind system was merely 1.02′′ (~ 3.6λ∕D). Since our observations operate at 10–12.5 μm, the diffraction-limited PSF core (0.28′′) is much wider than for high-contrast imaging at shorter wavelengths. At the small separations we are interested in, ADI-based reduction schemes are therefore impractical, since the self-subtraction imposed by such techniques become excessively large for modest amounts of field rotation. We thus adopted a recently introduced reduction technique called circularised PSF subtraction (Petit dit de la Roche et al. 2020), which is independent of field rotation and hence efficient for detecting companion signals at smaller angular separations. This technique involves creating a circularly symmetric stellar PSF by rotating the data frames through 360° in steps of 1° and taking the average of all such rotated versions of the frame. The stellar PSF thus created for each data frame is then subtracted from the original frame and the resulting frames are de-rotated to align the north upwards and east to the left. The final reduced image is then obtained by combining all the de-rotated reduced frames into a master image by means of a weighted sigma-clipped median function, choosing a threshold of 3σ, with the weights based on the standard deviation of each image. The circularised PSF subtraction technique is effectively equivalent to the classical approach of circular profile subtraction (Lafrenière et al. 2007), which has been commonly used in the past to reduce imaging data. In this reduction technique, the image is divided of into a sequence of narrow circles centred on the star. The mean (or median) of each circular area is then calculated and subtracted from all pixels within that area. We tested both approaches and found that they are indeed essentially equivalent. However, in addition to circular profile subtraction being computationally faster, in this circumstance, it also has a marginally better performance at the specific expected separation on ϵ Ind Ab, so we used it as the primary option for this study.

The result of the NEAR reduction is shown in Fig. 1b. The expected position of the companion with respect to the central star, as predicted from the best-fit ephemerides for the observation epoch, is also shown in the figure. As with NaCo data, no companion was directly visible in the final reduced image from the NEAR observation. The 5σ upper flux limits at different separations from the central star were calculated by fitting Gaussians at each such annulus with an FWHM of 1λ∕D and a peak equal to five times the standard deviation at that location. The photometric reference value used for this flux calculation was 5.68(± 5%) Jy, which is the flux in11.6 μm IRAS:12 band for ϵ Ind A as given in the VizieR photometry tool. The resulting sensitivity curve from the reduced image for NEAR observation is shown in Fig. 5, both from circularised PSF subtraction and the classical circular profile subtraction, with the latter being used for further analysis in this work. As can be predicted from the reduced image from Fig. 1b, negative residual images of the star from the nodding and chopping procedures worsen the contrast beyond ~10 au. However, we do not expect this to affect our results since the range of projected planet-star separation covered at NaCo and NEAR observation epochs from best-fit ephemerides is between 3.71 au (1.02′′) and 4.98 au (1.37′′).

As an alternative procedure, we also reduced the NEAR data using the TRAP algorithm, keeping in mind its usefulness at smaller separations relative to conventional ADI techniques. TRAP was run on a stacked dataset of all the frames from the three observation nights. The principal component fraction f for reduction was set to 0.3 just as in the case for NaCo reduction. However, the contrast obtained from TRAP reduction was less deep than that obtained via the classical circular profile subtraction technique in the range of separation we are interested in. Hence, for the work in this paper, we used the upper flux limits and contrast obtained via the latter technique for analysis. We accounted for attenuation in the obtained upper flux limits due to AGPM off-axis transmission in NEAR using the transmission curve provided in Maire et al. (2020a). The conversion of the resulting contrast to absolute magnitudes was calculated using the apparent magnitude of ϵ Ind A in the mid-infraredband WISE W3 (12.082 μm), m_W3 = 2.146, as a proxy to NEAR 10–12.5 μm band and the distance modulus for ϵ Ind A, μ = −2.1951.

Fig. 5 5σ upper flux limits from NEAR data reduction, obtained using the classical circular profile subtraction technique (blue solid line) and circularised PSF subtraction technique (grey dashed line). The flux limits from both techniques become worse beyond 10 au due to the stellar residuals from chopping and nodding features in the image.

4 Results and discussion

We derived mass detection limits for planetary companions as a function of distance from the central star from the obtained NaCo and NEAR magnitude limits using theoretical models for predicting the mass-luminosity relationships at different ages based on the brightness in each respective photometric band. For the purpose of our analysis, we used two different sets of models: the AMES-Cond (Allard et al. 2001; Baraffe et al. 2003) atmospheric and evolutionary model grids developed for T-dwarf and giant planet atmospheres, and the Morley et al. (2012, 2014) atmospheric models developed for T/Y dwarfs and Y dwarfs. AMES-Cond grids couple atmospheric models with evolutionary models to predict the evolution of effective temperature and luminosity as a function of age for different masses. The atmospheric model described by AMES-Cond grids assumes the dust grains to have settled down gravitationally below the photosphere, hence neglecting the dust opacity in the radiative transfer equation. These models are suitable for T_eff < 1400 K and are known to agree well with NIR photometry beyond 1 μm. We used the AMES-Cond grids to convert NaCo L′ absolute magnitudes to T_eff and mass (M_P) limits assuming ages of 1–5 Gyr. The temperature limits predicted by the AMES-Cond grids are very similar for all assumed ages from 1 to 5 Gyr; as an example, the predicted detectable T_eff for 5 Gyr as a function of increasing distance from the central star is shown in Figs. 6a,b corresponding to contrast curves at 50 and 90% TPR, respectively. In this and subsequent figures, we also show the minimum (1.02′′) and maximum (3.27′′) projected separation for the orbit of ϵ Ind Ab as determined by Feng et al. (2019), as well as the specific predicted separation at the epoch of the NaCo observations (1.37′′ ± 0.16′′, at a position of ⟨RA⟩ = 1.311′′ ± 0.179′′, ⟨Dec⟩ = −0.061′′ ± 0.392′′, relative to the star).

The corresponding mass detection limits at 50 and 90% TPR are shown in Figs. 7a,b, respectively, as a function of distance from the central star, for age assumptions of 1–5 Gyr. The companion mass predicted from Feng et al. (2019), ${3.25}_{-0.65}^{+0.39}$ M_J is also shown in the figures. Assuming an age of 1 Gyr for the system, at the predicted distance of 1.37′′ from the central star, the AMES-Cond grids predict a 5σ detection limit of T_eff = 207 K and M_P = 1.8 M_J (M${}_{L^{\prime }}=20.3$ mag) at 50% TPR and T_eff = 227 K and M_P = 2.1 M_J (M${}_{L^{\prime }}=19.7$ mag) at 90% TPR. At an age of 5 Gyr, the 5σ limits from the gridpredictions at 1.37′′ are T_eff = 208 K and M_P = 4.8 M_J at 50% TPR and T_eff = 234 K and M_P = 5.9 M_J at 90% TPR. For both 50 and 90% TPRs, the obtained mass limits from the AMES-Cond grids indicate that at an age of 1 Gyr, the companion should be detectable in NaCo L′ throughout its orbit. At 5 Gyr, ϵ Ind Ab does not fall within the mass detection limits in its range of projected orbital separation from the central star, at both 50 and 90% TPRs.

In atmospheric models like the one described by the AMES-Cond grids, the planet atmospheres are assumed to be free of dust opacity (and thereby clouds), the flux is predicted to arise in relatively deep layers of the atmosphere. By contrast, when cloud opacity is included in the radiative transfer equation, the depth into the atmosphere is restricted, thereby limiting the observed flux particularly at short NIR wavelengths. The Morley et al. (2012, 2014) atmospheric models include cloud opacity in the radiative-convective equilibrium model of brown dwarf atmospheres, with the Morley et al. (2012) grid suitable for T/Y dwarfs with surface temperatures between 400–1200 K and log(g) between 4 and 5.5; and theMorley et al. (2014) grid suitable for Y dwarfs with surface temperatures between 200–450 K and log(g) between 3 and 5. In addition to the cloud opacity for Na₂S, KCL, ZnS, MnS and Cr modelled in Morley et al. (2012), Morley et al. (2014) also includes H₂O and NH₃ that become significant condensates in the atmospheres of brown dwarfs cooler than 500 K. To consider the effect of clouds on the predicted mass limits, we used the Morley et al. (2012) atmospheric model for the background limited regime and Morley et al. (2014) atmospheric model for the contrast limited regime in our reduced NaCo data. Further, we constrained both these atmospheric models to solar metallicity and a log(g) = 4.0, which is the typical surface gravity of giant planets in the expected mass range as ϵ Ind Ab. The model is also described by a sedimentation efficiency parameter, f_sed, that decides the total amount of condensates assumed in each layer of atmosphere, directly affecting the predicted flux described by the model via Mie scattering. A higher f_sed factor describes optically thinner clouds in the atmosphere, that is, lesser vertical extension of the clouds and larger particle sizes. In this work, we assumed moderate optical thickness for the clouds and set the f_sed to a value of 5.0 in the model. Using the resulting constrained atmospheric model and based on the absolute NaCo L′ magnitudes obtained, we then predicted a T_eff for the companion at increasing orbital separations at both 50 and 90% TPRs, as shown in Figs. 6a,b.

In order to utilise the same evolutionary track as the AMES-Cond atmospheric and evolutionary model grids for the purpose of comparing the final mass limits, we used the predicted T_eff by Morley et al. (2012, 2014) atmospheric models to obtain a corresponding mass limit for the companion from the evolutionary grid underlying the AMES-Cond grids, assuming ages of 1–5 Gyr. The corresponding log(g) values obtained this way do not have any significant offset from the log(g) to which we constrained the Morley models. The mass detection limits thus obtained from contrast curves at 50 and 90% TPRs are shown in Figs. 7a,b, respectively. At 1 Gyr, at the predicted planet-star separation of 1.37′′, the Morley atmospheric models, coupled with evolutionary track from AMES-Cond grids, predict a 5σ detection limit of T_eff = 256 K and a mass M_P = 2.8 M_J at 50% TPR and T_eff = 273 K and a mass M_P = 3.25 M_J at 90% TPR for the planet. Assuming an age of 5 Gyr, at 1.37′′, the predicted T_eff and M_P are 256 K and 6.9 M_J at 50% TPR and 273 K and 7.8 M_J at 90% TPR. Based on the derived detection limits from Morley et al. (2012, 2014), ϵ Ind Ab should be detectable in NaCo L′ beyond 1.1′′ (4 au) at 50% TPR and beyond 1.37′′ (4.98 au) at90% TPR if the assumed age for the system is 1 Gyr. Similar to AMES-Cond grids predictions, the Morley grids also predict that the companion is undetectable in the NaCo L′ band within its projected separation range for an assumed age of 5 Gyr, for both 50 and 90% TPRs.

Table 1 summarises the 5σ T_eff, M_P detection limits for the companion obtained from NaCo L′ observations of ϵ Ind A for assumed age of 1 and 5 Gyr. The detectable mass limits from these observations are a significant improvement from those of previous imaging campaigns of ϵ Ind A; In comparison, at an age of 1 Gyr, using the Baraffe et al. (2003) models, Geißler et al. (2007) arrived at a detection limit of 21 M_J at separations of ≥1.3′′ and 16 ± 4 M_J at ≥ 3′′ from observations in the NaCo H and Ks bands, while Janson et al. (2009) arrived at a constraint of 5–20 M_J for the companion at separations of 10–20 au (~2.7′′–5.5′′) from observations in the NaCo NB4.05 narrow-band and L′ band.

We also applied the AMES-Cond atmospheric and evolutionary model grids and Morley et al. (2012, 2014) atmospheric models to NEAR data to derive detectable T_eff and mass limits for ϵ Ind Ab. Since the AMES-Cond atmospheric and isochronal grids currently available for VISIR photometric bands do not include the NEAR band, we calculated the brightness in this bandpass using the theoretical spectra for AMES-Cond available at the theoretical spectra web server developed by Spanish Virtual Observatory (SVO) (Bayo et al. 2008). We constrained the spectra to solar metallicities and a surface gravity of log(g) = 4.0. The NEAR band fluxes for ages 1–5 Gyr were obtained from the spectra by integrating it within the 10–12.5 μm bandpass and scaling the results by a factor of $\frac{R^2}{D^2}$. Here, D is the distance to ϵ Ind, and R is the planet radius corresponding to the spectral temperature obtained by interpolating the existing AMES-Cond grid for each respective age. To convert the obtained flux into magnitudes, we used a low-resolution mid-infrared spectrum of Vega from CASSIS (Lebouteiller et al. 2011) and integrated it over the NEAR band. We used this magnitude calibration and the distance modulus for ϵ Ind to calculate the absolute magnitudes in the NEAR band from the corresponding flux. We then interpolated the T_eff and NEAR magnitudes thus derivedfrom the AMES-Cond spectra into the grid points of the regular AMES-Cond grid. We used these new grid values to predict the detectable T_eff limits for ϵ Ind Ab in the NEAR band, with Fig. 8 showing the T_eff limits with increasing orbital separation for an assumed age of 5 Gyr.

The corresponding upper mass limits obtained for the companion in the NEAR band is shown in Fig. 9 as a function of orbital separation, for assumed ages of 1–5 Gyr. For an age of 1 Gyr, the AMES-Cond atmospheric and evolutionary model grids predict a 5σ detection limit of T_eff = 325 K and M_P = 4.8 M_J at the predicted planet-star separation of 1.02′′ in the NEAR epoch and for 5 Gyr, the 5σ limits obtained from the model at the same separation are T_eff = 338 K and M_P = 11.4 M_J. According to the predictions from AMES-Cond grids, the companion does not fall within the detectable range of T_eff and mass in NEAR band observations.

Similar to the analysis of NaCo data, we also applied Morley et al. (2012, 2014) atmospheric models to NEAR data to see how the limits change when clouds are included in the atmosphere. The current Morley grids do not include NEAR band magnitudes, and since the corresponding spectral files were unavailable, we chose the grid for WISE W3 band (12 μm) as a proxy for the NEAR band. Assuming ages of 1–5 Gyr for the system, we predicted the detectable T_eff limits from the obtained NEAR band magnitudes using the corresponding Morley grid. Figure 8 shows the 5σ T_eff limits for the planet with increasing orbital distance thus obtained from the Morley atmospheric model for an assumed age of 5 Gyr. As with NaCo data, we then used the evolutionary track for NEAR band AMES-Cond atmosphericand evolutionary model grids for the respective age to predict the final mass limits from the obtained T_eff limits. The derived upper mass limits for the planetary companion in the NEAR band using the Morley atmospheric model, coupled with evolutionary track from AMES-Cond grids, for different age assumptions is shown in Fig. 9. At 1 Gyr, the Morley grids predict a 5σ detection limit of T_eff = 284 K and M_P = 3.6 M_J at a separation of 1.02′′ from the star and at 5 Gyr, the predicted limits at the same separation are T_eff = 284 K and M_P = 8.4 M_J. In the same way as inferred from AMES-Cond limits, the Morley grids also predict that the expected mass of the companion, 3.25 M_J, is below the detectable mass range for NEAR observations. Table 2 summarises the 5σ T_eff, M_P detection limits for the companion obtained from NEAR 10–12.5 μm observations of ϵ Ind A for assumed age of 1 and 5 Gyr.

An interesting point to note from the above results is how in NaCo band predictions from the atmospheric models in AMES-Cond grids for T_eff and M_P are lower, hence better in this context, than the respective predictions from Morley et al. (2012, 2014) atmosphericmodels for the same age, while in the NEAR band Morley grids predict lower T_eff and M_P than AMES-Cond grids. This is a consequence of the Morley et al. (2012, 2014) atmospheric models being a bit redder at relevant temperatures, hence predicting higher flux at longer wavelengths than atmospheric models underlying AMES-Cond grids for the same temperature. This underlines the significance of using different models to predict the T_eff and mass limits of the companion to account for such model uncertainties while interpreting the results. Furthermore, it underlines the value of observing in multiple carefully selected wavelength bands as in this study since this provides redundancy against the existing model uncertainties and thereby yields much more robust detections or detection limits.

We note that we also attempted to use the Exo-REM atmospheric model (Baudino et al. 2015, 2017) as a third alternative model interpretation for placing constraints on the T_eff, and hence mass, of ϵ Ind Ab. The original Exo-REM model is essentially non-cloudy and valid for EGPs with T_eff = 500–2000 K. An upgraded model for Exo-REM, Charnay et al. (2018), includes both absorption and scattering by clouds in the atmosphere and is valid for EGPs and brown dwarfs with T_eff = 300–1800 K. However, as it turned out during our analysis, the detectable T_eff drops below 300 K for much of the separation range of ϵ Ind Ab. Since the Exo-REM model grids do not extend to such low temperatures, they are not applicable for the purpose of this analysis. Hence, we did not include them, but we do note that in the temperature range over which they overlap with the other tested model grids, both the cloudy and cloud-free Exo-REM models show good consistency with the atmospheric models underlying the AMES-Cond grids.

Fig. 6 Predicted 5σ temperature limits for ϵ Ind Ab from NaCo using contrast curves derived at 50 (panel a) and 90% TPR (panel b). The limits derived from the atmospheric models in AMES-Cond grids is shown as solid purple line and the limits from the Morley et al. (2012, 2014) atmospheric models is shown as solid green line. The age assumed for the above limits is 5 Gyr. The dashed black lines represent the minimum, maximum projected separation of the ϵ Ind Ab as well as the specific projected separation at the epoch of observation.

Fig. 7 5σ upper mass limits for ϵ Ind Ab derived from observations with NaCo L′, reduced at a TPR of 50 (panel a) and 90% (panel b), for different age assumptions. The solid lines represent limits obtained from AMES-Cond atmospheric and evolutionary grids and the dashed lines represent the limits from Morley et al. (2012, 2014) atmospheric models, coupled with the evolutionary track from AMES-Cond grids. The solid blue line represents the companion mass predicted from Feng et al. (2019) at an inclination of ${64.25}^{°}{}_{-6.09}^{+13.8}$, with the associated region of uncertainty in mass shown as the filled region between the blue dashed lines. The dashed black lines have the same meaning as in the previous figures.

Fig. 8 Predicted 5σ temperaturelimits for NEAR magnitudes as a function of increasing distance from the central star. The limits derived from AMES-Cond atmospheric and evolutionary grids is shown as solid purple line and the limits from Morley et al. (2012, 2014) atmospheric models is shown as solid green line. The age assumed for the above limits is 5 Gyr. The dashed black lines have the same meaning as in the previous figures.

Fig. 9 5σ upper mass limits for ϵ Ind Ab derived from NEAR band observations for different age assumptions. The symbols and colours have the same meaning as in the previous figures.

5 Conclusions

In this work, we present observations of ϵ Ind A using deep AO imaging at the VLT with NaCo in the L′ -band and NEAR in a dedicated 10–12.5 μm band, in order to attempt detection of the planetary companion ϵ Ind Ab. Assisted by the well-constrained ephemerides and mass of the planet, we can derive stringent constraints on its brightness from the imaging. We arrive at unprecedented sensitivities close to the bright star, with detectable planet temperatures as low as 200–300 K depending on the choice of atmospheric model between the AMES-Cond and the (Morley et al. 2014) models. From the corresponding planet mass detection limits obtained from AMES-Cond atmospheric and evolutionary models, the 3.25 M_J planet ϵ Ind Ab would be marginally detectable in the NaCo L′ images if the (uncertain) age of the ϵ Ind system is as low as 3 Gyr at a TPR of 50% and 2 Gyr at a TPR of 90%. On the other hand, the planet mass detection limits as obtained from Morley et al. (2014) models at both 50 and 90% TPR suggest that the planet will be detectable if the age of the system was 1 Gyr. As can be inferred from the non-detection of the planet in the images at its predicted location of 1.37′′, at 50% TPR, the AMES-Cond model grids thus hint at an age of 4 Gyr or above and the Morley et al. (2014) model grids hint at an age of 2 Gyr or above for the planet. At 90% TPR, the detectable planet mass and T_eff limits are slightly higher, with AMES-Cond model grids indicating an age of 3 Gyr or above but Morley et al. (2014) model grids predicting an age of 2 Gyr or above, consistent with the former case. Assuming that these two atmospheric models cover a somewhat accurate depiction of the companion’s atmosphere, both these cases definitely put a constraint of 2 Gyr for the lower age limit of the planetary candidate. The non-detection of the planet in these observations and the corresponding age estimation in this work is compatible with the indication of an older age for the system in the literature. The NaCo detection limits are more constraining than the NEAR limits according to both model sets, partly due to the high thermal background in NEAR and partly due to the less favourable separation of the planet at the NEAR epoch. It is also to be noted that while the NaCo data are contrast-limited at the expected position of the companion from the best-fit ephemerides, the NEAR data are essentially background limited as can be seen from the absence of any visible speckle structure around the expected planet position in Fig. 1b. A longer integration with NEAR may, hence, bring the contrast in this band further down, but for NaCo L′ it would not result in any significant improvement. However, while the difference in mass detection limits between NaCo and NEAR is very distinct in the AMES-Cond interpretation, it is much subtler in the Morley et al. (2014) interpretation, which predicts a worse sensitivity relative to AMES-Cond for NaCo, but a better relative sensitivity for NEAR. This underlines the important fact that the sensitivity in mass depends not only on the intrinsic instrumental sensitivity in different bands, but also on the model-predicted flux distribution among those bands. In this regard, having measurements in more than one high-sensitivity band as in our study adds considerable robustness against atmospheric uncertainties in the interpretation of the physical detection limits.

For securing an imaging detection of the planet at any realistic age, only a rather modest improvement in sensitivity (e.g. ~1 mag or less in L′-band) is required, particularly since the projected separation of the planet is expected to be more favourable for detection over the next few decades. The upcoming enhanced resolution imager and spectrograph (ERIS) (Kenworthy et al. 2018) instrument for the VLT is foreseen to be able to deliver the required sensitivity in the ~4 μm range. The mid-infrared instrument (MIRI) (Bouchet et al. 2015) on board the James Webb Space Telescope (JWST), which is currently set to launch in late 2021, offers broadband imaging with coronagraphy in a broad wavelength range of interest (10.6, 11.4, 15.5, 23 μm filters) at unmatched sensitivity (~1.5 μJy limiting sensitivity at 10σ for a 10 000s exposure at 11.3 μm, with a PSF FWHM of 0.36′′; Bouchet et al. 2015), aiding a possible detection of ϵ Ind Ab in MIR. The mid-infrared ELT imager and spectrograph (METIS) (Brandl et al. 2014, 2018; Quanz et al. 2015), an instrument on the 39 m European extremely large telescope (E-ELT) that is set to see its first light in the late 2020s, will offer mid-infrared imaging with AO and coronagraphy with a very high sensitivity (~0.1mJy at 10σ for a 1 h exposure at ~11 μm; Brandl et al. 2014) and contrast, enabling high-resolution spectroscopy for detailed atmospheric characterisation. In addition, further RV monitoring in combination with new astrometric data from future Gaia releases can place tighter constraints on the orbital parameters for ϵ Ind Ab, further aiding high-precision characterisation of a wide range of atmospheric and physical properties for this cold and very nearby giant planet.

Acknowledgements

M.J. gratefully acknowledges funding from the Knut and Alice Wallenberg Foundation. We made use of CDS and NASA/ADS services for the purpose of this study. Part of this work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. S.P.Q. and A.B. acknowledge the financial support from the SNSF. Part of this research received funding from the Fonds de la Recherche Scientifique – FNRS under Grant No F.4504.18, and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 819155). H.R.A.J. acknowledges support from the UK’s Science and Technology Facilities Council ST/R000905/1.

References

All Tables

All Figures

Fig. 1 Final detection probability map from NaCo (panel a) and the final reduced image from NEAR observations (panel b) of ϵ Ind A. The predicted position of the companion at the time of observation, along with the 1σ uncertainty in position, as inferred from the best-fit ephemerides, based on RV and astrometry data (Feng et al. 2019), is shown asan ellipse. The NaCo image was reduced using a combination of RSM and modified LOCI subtraction applied at suitable separations while the NEAR image was reduced using the circular profile subtraction technique. No companions were detected in either image.

In the text

	Fig. 2 Detection maps from NaCo reduction using RSM at 50 (panel a) and 90% TPR (panel b). Out of the 11 injected fake companions along the shown axis, four have been retrieved for the 50% case and nine for the 90% case. The detected signals correspond to a threshold of 0.0025, which is equal to the brightest speckle in the empty probability map.
In the text

Fig. 3 Contrast curves resulting from the different reduction algorithms applied to NaCo L′ observationsin this work. The y-axis shows the 5σ contrast limits obtained from TRAP, the normal LOCI, and modified LOCI subtraction techniques (shown as blue, green, and red dotted line), and the equivalent flux ratio or sensitivity limits obtained using RSM algorithm at 50 and 90% TPR as yellow and black dotted lines, respectively.

In the text

Fig. 4 Composite contrast curve from NaCo L′ observations, obtained from the combination of RSM and the modified LOCI subtraction techniques. The y-axis shows the 5σ contrast limits obtained from modified LOCI subtraction techniques (shown in red solid dots), beyond ~11 au and the equivalent flux ratio or sensitivity limits obtained using RSM algorithm corresponding to a TPR of 50% (blue solid dots) and 90% (green solid dots).

In the text

	Fig. 5 5σ upper flux limits from NEAR data reduction, obtained using the classical circular profile subtraction technique (blue solid line) and circularised PSF subtraction technique (grey dashed line). The flux limits from both techniques become worse beyond 10 au due to the stellar residuals from chopping and nodding features in the image.
In the text

Fig. 6 Predicted 5σ temperature limits for ϵ Ind Ab from NaCo using contrast curves derived at 50 (panel a) and 90% TPR (panel b). The limits derived from the atmospheric models in AMES-Cond grids is shown as solid purple line and the limits from the Morley et al. (2012, 2014) atmospheric models is shown as solid green line. The age assumed for the above limits is 5 Gyr. The dashed black lines represent the minimum, maximum projected separation of the ϵ Ind Ab as well as the specific projected separation at the epoch of observation.

In the text

Fig. 7 5σ upper mass limits for ϵ Ind Ab derived from observations with NaCo L′, reduced at a TPR of 50 (panel a) and 90% (panel b), for different age assumptions. The solid lines represent limits obtained from AMES-Cond atmospheric and evolutionary grids and the dashed lines represent the limits from Morley et al. (2012, 2014) atmospheric models, coupled with the evolutionary track from AMES-Cond grids. The solid blue line represents the companion mass predicted from Feng et al. (2019) at an inclination of ${64.25}^{°}{}_{-6.09}^{+13.8}$, with the associated region of uncertainty in mass shown as the filled region between the blue dashed lines. The dashed black lines have the same meaning as in the previous figures.

In the text

Fig. 8 Predicted 5σ temperaturelimits for NEAR magnitudes as a function of increasing distance from the central star. The limits derived from AMES-Cond atmospheric and evolutionary grids is shown as solid purple line and the limits from Morley et al. (2012, 2014) atmospheric models is shown as solid green line. The age assumed for the above limits is 5 Gyr. The dashed black lines have the same meaning as in the previous figures.

In the text

	Fig. 9 5σ upper mass limits for ϵ Ind Ab derived from NEAR band observations for different age assumptions. The symbols and colours have the same meaning as in the previous figures.
In the text