A faint companion around CrA-9: protoplanet or obscured binary?
V. Christiaens, M.-G. Ubeira-Gabellini, H. Cánovas, P. Delorme, B. Pairet, O. Absil, S. Casassus, J. H. Girard, A. Zurlo, Y. Aoyama, G-D. Marleau, L. Spina, N. van der Marel, L. Cieza, G. Lodato, S. Pérez, C. Pinte, D. J. Price, M. Reggiani
MMNRAS , 1–22 (2021) Preprint 23 February 2021 Compiled using MNRAS L A TEX style file v3.0
A faint companion around CrA-9: protoplanet or obscured binary?
V. Christiaens ★ , M.-G. Ubeira-Gabellini , H. Cánovas , P. Delorme , B. Pairet ,O. Absil , S. Casassus , J. H. Girard , A. Zurlo , , Y. Aoyama , , G-D. Marleau , , ,L. Spina , N. van der Marel , L. Cieza , , G. Lodato , S. Pérez , C. Pinte , ,D. J. Price , M. Reggiani School of Physics and Astronomy, Monash University, Clayton Vic 3800, Australia Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano MI, Italy Aurora Technology for ESA/ESAC, Camino bajo del Castillo s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain Institut de Planétologie et d’Astrophysique de Grenoble, Université Grenoble Alpes, 38058 Grenoble, France ISPGroup, ELEN/ICTEAM, UCLouvain, Belgium Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile Space sciences, Technologies & Astrophysics Research (STAR) Institute, Université de Liège, Allée du Six Août 19c, B-4000 Sart Tilman, Belgium Space Telescope Science Institute, 3700 San Martin Dr. Baltimore, MD 21218, USA Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile Escuela de Ingeniería Industrial, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile Institute for Advanced Study, Tsinghua University, Beijing 100084, People’s Republic of China Department of Astronomy, Tsinghua University, Beijing 100084, People’s Republic of China Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germany Physikalisches Institut, Universität Bern, Gesellschaftsstr. 6, CH-3012 Bern, Switzerland Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany INAF-Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, Padova, Italy, 35122-I Physics & Astronomy Department, University of Victoria, 3800 Finnerty Road, Victoria, BC, V8P 5C2, Canada Departamento de Física, Universidad de Santiago de Chile. Avenida Ecuador 3493, Estación Central, Santiago, Chile Institute of astrophysics, KU Leuven, Celestijnlaan 200D, 3001 Leuven, Belgium
Accepted 2021 February 16. Received 2021 February 15; in original form 2020 December 01
ABSTRACT
Understanding how giant planets form requires observational input from directly imaged pro-toplanets. We used VLT/NACO and VLT/SPHERE to search for companions in the transitiondisc of 2MASS J19005804-3645048 (hereafter CrA-9), an accreting M0.75 dwarf with anestimated age of 1–2 Myr. We found a faint point source at ∼ (cid:48)(cid:48) separation from CrA-9 ( ∼ 𝜎 significance. The near-IR absolute magnitudes of the object point towards a planetary-masscompanion. However, our analysis of the 1.0–3.8 𝜇 m spectrum extracted for the companionsuggests it is a young M5.5 dwarf, based on both the 1.13- 𝜇 m Na index and comparison withtemplates of the Montreal Spectral Library. The observed spectrum is best reproduced with higheffective temperature (3057 + − K) BT-DUSTY and BT-SETTL models, but the correspondingphotometric radius required to match the measured flux is only 0 . + . − . Jovian radius. Wediscuss possible explanations to reconcile our measurements, including an M-dwarf compan-ion obscured by an edge-on circum-secondary disc or the shock-heated part of the photosphereof an accreting protoplanet. Follow-up observations covering a larger wavelength range and/orat finer spectral resolution are required to discriminate these two scenarios.
Key words: protoplanetary discs – planet-disc interactions – techniques: image processing –planets and satellites: formation ★ E-mail: [email protected]
The classical debate on the formation of giant planets confronts core accretion (Mizuno 1980; Pollack et al. 1996) to gravitationalinstability (Boss 1998; Kratter & Lodato 2016). While the majority © a r X i v : . [ a s t r o - ph . E P ] F e b V. Christiaens et al. of the detected population of short-orbit mature exoplanets appearsconsistent with predictions from core accretion models (e.g. Winn& Fabrycky 2015; Mordasini 2018), it is unclear whether the prop-erties of young giant planets that have been directly imaged at largeorbital separations are also consistent with formation through coreaccretion (e.g. the HR 8799 planets, 𝛽 Pic b and c, HIP 65426 b;Marois et al. 2008; Lagrange et al. 2009; Bonnefoy et al. 2013;Marleau & Cumming 2014; Chauvin et al. 2017; Marleau et al.2019a; Nowak et al. 2020). How did these adolescent 5–12 Jupitermass ( 𝑀 𝐽 ) planets found at up to ∼
100 au separation form in thefirst place? If similar planetary-mass companions are also foundat large separations at very young ages ( ∼ transition discs , constitute primetargets to search for nascent giant planets, since they may be carvingthe cavity (e.g. Espaillat et al. 2014; Casassus 2016; Owen 2016;van der Marel et al. 2021).High-contrast imaging in IR is one of the most powerful tech-nique to detect those young companions (e.g. Absil & Mawet 2010;Bowler 2016). A particularly suited observing strategy to reachhigh contrast is angular differential imaging (ADI; Marois et al.2006). When coupled with an appropriate post-processing algo-rithm such as principal component analysis (PCA; Amara & Quanz2012; Soummer et al. 2012), this technique can efficiently modeland suppress the bright stellar halo of the star, while preserving thatof the planet. Nevertheless, in the presence of a bright circumstellardisc, aggressive ADI filtering can create point-like artefacts whichcan be confused with substellar companions (e.g. Milli et al. 2012;Rich et al. 2019; Currie et al. 2019). This has led to a number ofprotoplanet detection claims whose authenticity has subsequentlybeen debated in the recent years (e.g Quanz et al. 2013; Sallum et al.2015; Guidi et al. 2018). To add to the confusion, even faint com-panions imaged at a location external to the circumstellar disc canalso be misclassified. The IR magnitudes of companions FW Tau Cand CS Cha B suggested a planetary mass (Kraus et al. 2014; Ginskiet al. 2018), however recent studies have shown that they were morelikely to be obscured M-dwarf companions (Wu & Sheehan 2017;Haffert et al. 2020).So far, the only confirmed detection of protoplanets was madein the transition disc of PDS 70, with multiple independent detec-tions in the IR (Keppler et al. 2018; Müller et al. 2018; Christiaenset al. 2019a,b), in the H 𝛼 line (Haffert et al. 2019), and at sub-mmwavelengths (Isella et al. 2019). This further motivates the searchfor young companions in other discs harbouring large cavities. Inparticular, a statistically significant number of detections at veryyoung age could constrain the planet formation mechanisms thatare indeed at work, and the connection with their natal disc.In this work we focus on 2MASS J19005804-3645048 (here-after CrA-9 , as in Peterson et al. 2011; Romero et al. 2012; Caz-zoletti et al. 2019), a young accreting T-Tauri star in the CoronaAustralis (CrA) molecular cloud, surrounded by a transition disc.We report the discovery and characterisation of a faint point sourceat 0 . (cid:48)(cid:48) Table 1.
Characteristics of CrA-9 and its protoplanetary disc.Parameter Value Reference2MASS name J19005804-3645048Right ascension 19 h m s .
044 1Declination − ◦ (cid:48) . (cid:48)(cid:48)
883 1Distance [pc] 153 . ± . ± 𝑇 eff [K] 3720 ± ( a ) 𝑔 ) 3 . ± . ( b ) Luminosity 𝐿 ★ [ 𝐿 (cid:12) ] 0.46 2Age [Myr] 1–2 Sec. 2.2Mass [ 𝑀 (cid:12) ] 0 .
45 3Accretion rate [Log( 𝑀 (cid:12) yr − )]] -8.6 2Li EW [Å] 0.48 2 𝐴 V [mag] 1.8–2.1 4, 5Disc luminosity [Log( 𝐿 d / 𝐿 ★ )] -2.4 2Dust mass [ 𝑀 ⊕ ] 3 . ± .
12 3
Notes : ( a ) Based on the empirical relation to convert from spectral typeto effective temperature in Herczeg & Hillenbrand (2014). ( b ) Based oneffective temperature and age estimates, and the isochrones of either Tognelliet al. (2011) or Baraffe et al. (2015).
References : (1) Gaia Collaboration et al. (2018) (2) Romero et al. (2012);(3) Cazzoletti et al. (2019); (4) Dunham et al. (2015); (5) van der Marel et al.(2016). the point source. In Section 5, we discuss its possible nature. Finally,we summarise our conclusions in Section 6.
Table 1 summarises the known physical properties of CrA-9. CrA-9is located at the edge of the R CrA dark cloud (also referred to as the
Coronet ; Taylor & Storey 1984), a highly obscured and very youngregion of the CrA molecular cloud (Figure 1a and b; Gutermuthet al. 2009; Peterson et al. 2011; Bresnahan et al. 2018). The GaiaDR2 parallax for CrA-9 corresponds to a distance of 153 . ± . ± ± molecular band (Cruz &Reid 2002). Using the empirical relation in Herczeg & Hillenbrand(2014), this converts to an effective temperature of 3720 ± 𝐿 (cid:12) .An extinction of 𝐴 𝑉 ∼ . 𝐴 𝑉 ∼ . − for the H 𝛼 line, which convertsinto an accretion rate (cid:164) 𝑀 ≈ . × − 𝑀 (cid:12) yr − (Natta et al. 2004). MNRAS000
Coronet ; Taylor & Storey 1984), a highly obscured and very youngregion of the CrA molecular cloud (Figure 1a and b; Gutermuthet al. 2009; Peterson et al. 2011; Bresnahan et al. 2018). The GaiaDR2 parallax for CrA-9 corresponds to a distance of 153 . ± . ± ± molecular band (Cruz &Reid 2002). Using the empirical relation in Herczeg & Hillenbrand(2014), this converts to an effective temperature of 3720 ± 𝐿 (cid:12) .An extinction of 𝐴 𝑉 ∼ . 𝐴 𝑉 ∼ . − for the H 𝛼 line, which convertsinto an accretion rate (cid:164) 𝑀 ≈ . × − 𝑀 (cid:12) yr − (Natta et al. 2004). MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? This corresponds to an accretion luminosity 𝐿 acc ≈ . 𝐿 (cid:12) , henceabout 8% of the stellar luminosity. Romero et al. (2012) measured an equivalent width (EW) of 0.48Å for the Li I line (6707 Å). Considering the effective temperature ofthe star, the non-depletion of lithium points towards an age youngerthan ∼ 𝜌 Ophiucus dark cloud. Both the measured 0.48 Å EW and theLi I EWs of other members of the CrA cloud (e.g. Sicilia-Aguilaret al. 2008) are consistent to what is seen in other stellar formingregions with ages ∼ (cid:46) ±
150 K) andluminosity ( ∼ . 𝐿 (cid:12) ) inferred for CrA-9 assuming a distanceof 150 pc (Romero et al. 2012), hence consistent with all otherestimates discussed above. The SED of CrA-9 is characteristic of a protoplanetary disc with alarge dust cavity (Romero et al. 2012; van der Marel et al. 2016), i.e.a transition disc . Based on both the positive slope of the IR excessand the high stellar accretion rate, Romero et al. (2012) recognisedit as one of the transition discs of their sample whose properties arethe most compatible with giant planet(s) dynamically carving thedust depletion. van der Marel et al. (2016) inferred a cavity radiusof ∼
14 au based on SED modeling.Sicilia-Aguilar et al. (2008) had already noted the larger frac-tion of transition discs in the R CrA cloud core ( ∼ 𝐿 d / 𝐿 ★ ≈ − . ;Romero et al. 2012) also suggests an evolutionary stage earlier thandebris discs (typically 𝐿 d / 𝐿 ★ < − ; e.g. Cieza et al. 2010). Thisis also consistent with the fact that the star is still accreting. The disc has been detected by ALMA 1.3mm continuum ob-servations, as part of a survey of the whole CrA cloud, with anestimated dust disc radius of ∼ . (cid:48)(cid:48)
39 (Cazzoletti et al. 2019). Con-sidering standard assumptions on dust opacity at mm wavelengths(Beckwith et al. 1990) and a constant dust temperature of 20 K, themeasured 1.3-mm flux of 5 . ± .
16 mJy translates to 3 . ± . 𝑀 ⊕ of dust for an optically thin disc (Cazzoletti et al. 2019). Under thestandard (but highly uncertain) assumption of a gas-to-dust massratio of 100:1 (Williams & Cieza 2011), this corresponds to a discmass of ∼ . 𝑀 𝐽 . We observed CrA-9 with VLT instruments NACO (Rousset et al.2003; Lenzen et al. 2003) and SPHERE (Beuzit et al. 2008; Claudiet al. 2008; Dohlen et al. 2008) at four different epochs (ESO pro-grams 099.C-0883, 0101.C-0924, 0103.C-0865, and 103.2036.001,respectively). Owing to the faintness of the source, all observa-tions were obtained without coronagraph. Table 2 summarises theobservations. 𝐻 -band polarimetric dataset We observed CrA-9 with the NACO instrument (Rousset et al. 2003;Lenzen et al. 2003) at the VLT on 8 June 2019 in service mode(ESO programme 0103.C-0865). The observations were performedin polarimetric mode using the broad-band NACO 𝐻 -band filter( 𝜆 𝑐 = . 𝜇 m). In this observing mode a half-wave plate (HWP)rotates the polarisation plane of the incoming light before a Wol-laston prism splits the light into two orthogonally polarised beamsthat are projected on different detector regions. The CCD pixel sizewas set to 0 . (cid:48)(cid:48)
027 px − , the readout mode to Double RdRstRd , andthe detector mode to
HighDynamic , while we used the K dichroicthat splits the incoming light between the CONICA system and thewavefront-sensor.The observations consisted in multiple polarimetric cycleswhere each cycle contains four datacubes, one per HWP position an-gle (at 0 ◦ , . ◦ , ◦ , and 67 . ◦ , measured on sky east from north).We used detector integration times (DIT’s) of 0 . . . . . (cid:48)(cid:48) ± . (cid:48)(cid:48)
13. Standard cali-brations including darks and flat fields were provided by the ESOobservatory.The two simultaneous, orthogonally polarised images recordedon the detector when the HWP is at 0 ◦ ( ◦ ) were subtracted toproduce the Stokes parameter 𝑄 + ( 𝑄 − ) . This process was repeatedfor the 22 . ◦ ( . ◦ ) angles to produce the Stokes 𝑈 + ( 𝑈 − ) images.The total intensity (Stokes I) was computed by combining all theimages. We used the imaging polarimetry pipeline described byCanovas et al. (2011) and Canovas et al. (2015) to process the rawdata. Each science frame was dark-current subtracted and flat-fieldcorrected. Hot and dead pixels were identified with a 𝜎 -clippingalgorithm and masked out using the average of their surroundinggood pixels. The two images recorded in each science frame werealigned with an accuracy of 0.05 pixels. This process was applied toevery science frame resulting in a datacube for each Stokes 𝑄 ± , 𝑈 ± parameter. We applied the double-difference method as describedin Canovas et al. (2011) to correct for instrumental polarisation the MNRAS , 1–22 (2021)
V. Christiaens et al.
Table 2.
Summary of the observations on CrA-9 used in this work.Date Program Instrument Filter Mode Plate scale DIT a NDIT b NEXP c T dint < 𝛽 > e Δ PA f Notes[mas px − ] [s] [min] [ (cid:48)(cid:48) ] [ ◦ ]2017-08-26 099.C-0883 NACO 𝐿 (cid:48) ADI 27.2 0.1 453 62 44.7 0.52 25.3 non-sat.2018-07-02 0101.C-0924 NACO 𝐿 (cid:48) ADI 27.2 0.1 453 124 88.0 0.95 42.3 non-sat.2019-07-02 0103.C-0865 NACO 𝐻 PDI 27.15 0.8 80 56 56 0.57 – non-sat.2019-09-28 103.2036.001 IFS
𝑌 𝐽 𝐻
ADI 7.46 16 11 32 80.3 0.75 30.5 non-sat.2019-09-28 103.2036.001 IRDIS 𝐾
12 ADI 12.256 8 22 32 71.7 0.75 30.5 sat. core2019-09-28 103.2036.001 IRDIS 𝐾
12 ADI 12.256 2 24 2 1.1 0.75 31.8 non-sat. a Detector integration time. b Number of co-add images in each exposure. c Number of exposures. d Total integration times (excluding overheads), calculated after bad frame removal as explained in Section 3.1.3. e Average seeing at 𝜆 = f Parallactic angle variation achieved during the observed sequence. final, median-combined images. Finally, we derived the polarisedintensity ( 𝑃 𝐼 = √︁ 𝑄 + 𝑈 ) and the 𝑄 𝜙 and 𝑈 𝜙 images (see Schmidet al. 2006). 𝐿 (cid:48) -band datasets Both the 2017 and 2018 NACO datasets were acquired using the 𝐿 (cid:48) filter ( 𝜆 ∼ . 𝜇 m). The first dataset was obtained in servicemode in excellent conditions on 26 August 2017 (stable DIMMseeing ∼ . (cid:48)(cid:48) ∼ . (cid:48)(cid:48) ∼ (cid:48)(cid:48) . To compensate, we acquired atwice longer integration than at the first epoch (88 min vs. 45 min).Both observations were obtained in cube mode, which allows oneto save each individual co-add image instead of median-combiningthem. We opted for a pupil-tracking observing strategy to enableangular differential imaging (ADI; Marois et al. 2006), achieving25 ◦ and 42 ◦ field rotation in the 2017 and 2018 datasets, respec-tively. We used a 2-point dithering pattern in the two good quad-rants of NACO’s detector, excluding the bottom-left and top-rightquadrants to avoid bad columns and higher dark current noise, re-spectively. The stellar PSF did not saturate the detector in eitherobservation.We reduced both datasets in the same way, using a custom-made pipeline built on Python routines from the Vortex ImageProcessing package (vip ; Gomez Gonzalez et al. 2017). In brief,we (i) found the approximate stellar position in each science cubeand record the quadrant; (ii) subtracted an estimate of the sky back-ground for each image using images where the star is in a differentquadrant; (iii) flat-fielded each image; (iv) corrected for NaN val-ues and bad pixels; (v) found the stellar centroid with a 2D Moffatfit and shift all images to place the star on the central pixel; (vi)combined all images in a master cube and compute the associatedparallactic angles; (vii) identified and remove bad frames; (viii)median-combined together sets of 10 (resp. 16) consecutive imagesin the 2017 (resp. 2018) dataset; (ix) measured the FWHM in themedian image and the stellar flux in each image of the cube; and(x) finally post-processed the calibrated cube using median-ADI(Marois et al. 2006) and principal component analysis (PCA-ADI;Amara & Quanz 2012; Soummer et al. 2012). Appendix A givesthe details of each step. Available from https://github.com/vortex-exoplanet/vip . We searched for NACO 𝐿 (cid:48) observations of standard stars inthe ESO archive in order to provide an absolute photometric cali-bration of CrA-9 in the 𝐿 (cid:48) band. No standard star was observed thesame night as our observations of CrA-9, therefore we extended oursearch to standard stars observed in the same airmass and seeingconditions. We considered standard stars HD 205772 and HD 75223observed on 27 August 2017 and 4 December 2018, for the NACO2017 and 2018 CrA-9 datasets, respectively. We applied steps (i) to(v) of our NACO reduction pipeline for sky+dark subtraction, flat-fielding, bad pixel correction and recentering of the standard starsdata (Appendix A). Visual inspection of the PSF of the STD starssuggests similar Strehl ratios as achieved for CrA-9, for each respec-tive pair of observations. We measured the average flux (in ADUs)in a 1-FWHM aperture, and used the physical 𝐿 (cid:48) fluxes tabulatedin van der Bliek et al. (1996) to compute zero points. The absolute 𝐿 (cid:48) fluxes calculated for CrA-9 in the 2017 and 2018 datasets are1 . ± . × − W m − 𝜇 m − and 1 . ± . × − W m − 𝜇 m − respectively, where we considered a 10% relative uncertaintybased on our procedure (e.g. possible small differences in achievedStrehl ratios for the standard stars and CrA-9 observations). Despitethe uncertainties involved with our procedure, we note that thesedifferent absolute fluxes for the star appear to compensate exactlythe discrepant contrasts measured at the two NACO epochs for thefaint point source (i.e. they lead to the same 𝐿 (cid:48) flux Section 4.2.2). We observed CrA-9 on 28 September 2019 with VLT/SPHERE(Beuzit et al. 2008) in IRDIFS-EXT mode, i.e. with both the IFS andIRDIS sub-instruments acquiring images simultaneously (Claudiet al. 2008; Dohlen et al. 2008). IFS has a spectral resolution of 54,covering the
𝑌 𝐽𝐻 bands. It acquired 32 datacubes with 11 co-adds,each containing 39 spectral channels ranging from 0.95 to 1.68 𝜇 m.The stellar PSF did not saturate in any of the spectral channels.The same number of datacubes was obtained with IRDIS in the 𝐾 𝐾 𝜆 ≈ . 𝜇 m and 2 . 𝜇 m, respectively), withan integration time of 8s. The core of the stellar PSF saturated inthat sequence. We also acquired two datacubes at the beginning andend of the observation with the integration time set to 2s in orderto measure the unsaturated stellar flux. Conditions were averageand variable throughout the sequence (average seeing of ∼ . (cid:48)(cid:48) 𝑅 ≈ . MNRAS000
𝑌 𝐽𝐻 bands. It acquired 32 datacubes with 11 co-adds,each containing 39 spectral channels ranging from 0.95 to 1.68 𝜇 m.The stellar PSF did not saturate in any of the spectral channels.The same number of datacubes was obtained with IRDIS in the 𝐾 𝐾 𝜆 ≈ . 𝜇 m and 2 . 𝜇 m, respectively), withan integration time of 8s. The core of the stellar PSF saturated inthat sequence. We also acquired two datacubes at the beginning andend of the observation with the integration time set to 2s in orderto measure the unsaturated stellar flux. Conditions were averageand variable throughout the sequence (average seeing of ∼ . (cid:48)(cid:48) 𝑅 ≈ . MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? up to a factor ∼ ∼ ◦ field rotation was achieved throughoutthe pupil-stabilised sequence.We implemented two new pipelines to reduce our non-coronagraphic IFS and IRDIS data, respectively. Both of them aredivided in three parts: basic calibration, pre-processing and post-processing. For both pipelines, the basic calibration relied mostly onESO’s Common Pipeline Library esorex recipes (version 3.13.2),while both the pre- and post-processing parts made use of routinesfrom the vip package.Our reduction pipeline for IRDIS data consists of: (i) sky back-ground subtraction; (ii) flat-fielding; (iii) bad pixel correction; (iv)centering based on 2D Moffat fits of the stellar halo; (v) trimmingof bad frames; (vi) correction for the anamorphism present in theIRDIS images (Maire et al. 2016); and (vii) median-ADI, PCA-ADIand sPCA post-processing (Absil et al. 2013). Appendix B detailsthe different steps of the pipeline.Our IFS reduction pipeline involves more steps than for IRDIS,owing to both the complexity of the IFS instrument and the iden-tification of some sub-optimal features in the esorex calibrationrecipes. In short, the pipeline first computes all master calibrationfiles (darks, coloured and white detector flat fields, spectra positions,total instrument flat, wavelength solution, IFU flat) and uses them toreduce the science images and convert them into spectral cubes of39 channels. Next, vip routines deal with the bad pixel correction,centering on the star in each frame, bad frame identification and re-moval, and anamorphism correction. Finally the pipeline also relieson vip for post-processing of the cubes leveraging on either spec-tral differential imaging (SDI; Sparks & Ford 2002) and/or angulardifferential imaging (ADI). More specifically, we used PCA-SDI onindividual spectral cubes, PCA-ADI on 3D cubes for each spectralchannel sampled in the temporal dimension, and PCA-ASDI on the4D cubes leveraging on both SDI and ADI with a single PCA li-brary. Details on each step of the pipeline and on the post-processingare provided in Appendix C. In particular, we detail the calibrationsteps that enabled us to mitigate bright stripes in the final images.We double-checked the performance of our IRDIS and IFSreduction pipelines by comparing our calibrated cubes to the outputsof an independent calibration made by the SPHERE data center(Delorme et al. 2017a; Galicher et al. 2018). We found consistentastrometric and contrast estimates for the point source in the IRDISdata (Section 4.2). However the SPHERE data center calibration ofthe IFS data led to bright stripes in the final post-processed images.We also identified sub-optimal steps regarding the dark subtractionand distortion correction, which are yet to be implemented. Wetherefore favoured the results obtained by our pipeline for the restof this work.Contrary to our NACO data, we do not expect to achieve a goodabsolute flux calibration of our SPHERE stellar flux measurementsbased on STD stars. This is because CrA-9 is significantly fainterthan the nominal R mag limit for the visible wavefront sensor ofSPHERE to provide a good AO correction (R ∼ > 𝑅 band,and are hence expected to achieve a better Strehl ratio. We thereforeconsidered literature flux measurements of CrA-9 itself for absoluteflux calibration (Section 4.3). negfc module of vipThe negative fake companion technique (NEGFC; e.g. Lagrangeet al. 2010; Marois et al. 2010) enables to extract reliable astrometryand photometry for faint point sources found in images obtained us-ing ADI-related post-processing algorithms. The principle of negfcis to inject a negative companion in the calibrated cube (i.e. beforeADI post-processing), and find the position and flux that minimisethe residuals in the final ADI image in an aperture centered on thelocation of the companion candidate. This forward modeling ap-proach allows us to alleviate the biases that would affect astrometricand photometric estimates made directly in the final ADI image, i.e.geometric biases and flux losses. In this work, we have updated thenegfc module implemented in vip (Wertz et al. 2017; Gomez Gon-zalez et al. 2017) in order to improve the astrometric and contrastestimates of the faint point source presented in Section 4.1.The negfc module in vip relies on PCA-ADI in a single 3-FWHM wide annulus including the companion candidate. By de-fault, the default PCA algorithm used in NEGFC does not considera threshold in PA to build the PCA library (for computation effi-ciency). We have now added the option to use a threshold in PA(which is used in this work for the non-saturated IRDIS dataset).Three consecutive steps are involved for refined estimates of thecompanion’s radial separation, PA and flux: (1) a grid search onthe negative flux alone, at companion coordinates provided by theuser; (2) a Nelder-Mead downhill simplex on the three free param-eters, using the estimates in step 1 (Nelder & Mead 1965); and(3) a Markov Chain Monte Carlo (MCMC) algorithm sampling theprobability distribution of the companion’s 3D parameter space,using the simplex result as initial guess. The MCMC algorithm re-lies on emcee (Foreman-Mackey et al. 2013), which is a Pythonimplementation of the affine-invariant ensemble sampler proposedin Goodman & Weare (2010), and allows to infer uncertainties oneach of the three parameters.We applied three changes to the MCMC algorithm comparedto the version presented in Wertz et al. (2017): (i) we now use a newexpression for the log-probability provided to the MCMC sampler;(ii) we now check the convergence of the MCMC algorithm basedon the integrated auto-correlation time instead of a Gelman-Rubintest (Gelman & Rubin 1992); and (iii) we have now added theoption to inject different negative companion fluxes in the differentframes of the ADI cube, according to weights reflecting varyingobserving conditions throughout the ADI sequence. These changesare detailed in Appendix D. specfit : A new module for the spectral characterisationof point sources We implemented specfit, a new vip (Gomez Gonzalez et al. 2017)module which provides the tools to perform the spectral characteri-sation of directly imaged companions in a Bayesian framework. Thecore routine of the module, mcmc_spec_sampling , is a wrapper ofthe emcee package (Foreman-Mackey et al. 2013), adapted to sam-ple the probability distribution of the free parameters associated tothe models that are fitted to the observed spectrum. Any grid ofmodel can be used for the fit as long as a snippet function to readinput grid files is provided as argument to mcmc_spec_sampling .Apart from the parameters associated to the model grid suchas the effective temperature ( 𝑇 𝑒 ) and surface gravity (log ( 𝑔 ) ), addi-tional free parameters include (i) the photometric radius ( 𝑅 ) usedto scale the model along with the provided distance to the sys-tem; (ii) optionally the optical extinction 𝐴 𝑉 , treated as in Cardelli MNRAS , 1–22 (2021)
V. Christiaens et al.
Herschel - 250 µm Herschel - 70 µmb)a) 1 pc 0.2 pc
Coronet a r c s ec
50 auPCA-ADI (npc=10)
SPHERE/IRDIS K1K2 (2019-09-30) a r c s ec
50 auPCA-SDI (npc=1)
SPHERE/IFS YJH (2019-09-30) a r c s ec
50 auPCA-ADI (npc=23)
NACO L' (2018-07-02) a r c s ec
50 aumedian-ADI
NACO L' (2017-08-26) c) d) e) f)spirals?
Figure 1. a–b)
Herschel images showing the CrA-9 system (circled) within the Corona Australis star forming region (Bresnahan et al. 2018). CrA-9 is locatedwithin 1 pc projected separation from the dark cloud (
Coronet ), hence suggesting a very young age. c–f)
Images of the CrA-9 system obtained with VLT/NACO(c and d) and VLT/SPHERE (e and f) after subtraction of stellar emission using either angular differential imaging (c, d and f) or spectral differential imaging(e). The color scale of panels c-f) is linear and the cuts used are min/max except for panel e), where we used a scale spanning 0.1–99.9 percentiles in order tohighlight a tentative spiral pattern possibly connected to the companion. et al. (1989); (iii) optionally the ratio of total to selective extinction 𝑅 𝑉 , set by default to the diffuse interstellar medium (ISM) value 𝑅 𝑉 = . L( 𝐷 | 𝑀 ) = − 𝜒 (1) 𝜒 = (cid:2) W ( F obs − F mod ) 𝑇 (cid:3) C − (cid:2) W 𝑇 ( F obs − F mod ) (cid:3) (2)where W is the vector of weights 𝑤 𝑖 ∝ Δ 𝜆 𝑖 / 𝜆 𝑖 , with Δ 𝜆 𝑖 the FWHMof the filter (for IRDIS and NACO points) or spectral channel width(for IFS points); F obs and F mod are the fluxes of the observed andmodel spectra; and C is the spectral covariance matrix. W is normalised so that (cid:205) 𝑖 𝑤 𝑖 = 𝑁 , where N is the number ofpoints in the spectrum. The inclusion of W in the expression of thelog-likelihood makes it different to that used in recent MCMC-basedSED modeling implementations (e.g. Wang et al. 2020; Stolker et al.2020a; Wang et al. 2021). The motivation behind the use of W is toassign a weight proportional to the amount of “spectral information”contained by each point. Without W , all points would contributein the same way to the likelihood. For our spectral sampling ofthe point source, this would bias the algorithm in trying to betterreproduce the 𝑌 𝐽𝐻 points, where a higher density of measurementsis available, at the expense of the 𝐾 𝐾 𝐿 (cid:48) photometric points,although the latter cover a larger bandpass.The model flux points 𝐹 mod , i were obtained after convolutionwith the filter of the respective instrument they are compared to.In the case of the IFS channels, we considered a 17.33-nm FWHMbased on the specifications of the IFS prism provided in the ESOmanual, while for the 𝐾 𝐾 𝐿 (cid:48) points we used the filtertransmission curves provided by the observatory. The values of 𝐶 𝑖 𝑗 for 𝑖 and 𝑗 <
39 are calculated as in Delorme et al. (2017b) on thePCA-ADI images obtained with the different IFS spectral channels,and 𝐶 𝑖 𝑗 = 𝛿 𝑖 𝑗 for 𝑖 or 𝑗 >
39 (i.e. the IRDIS and NACO points),where 𝛿 𝑖 𝑗 is the Kroenecker symbol.The mcmc_spec_sampling routine allows to infer the mostlikely parameter values for a given parametric model grid. How-ever, for fits to non-parametric libraries, we have implemented best_fit_tmp , a routine to search for the most similar templateto an input spectrum, which is agnostic of the chosen spectral li- MNRAS , 1–22 (2021) rA-9 b, B or BKG? arcsec a r c s e c Stokes P arcsec Q arcsec U Figure 2.
Stokes 𝑃 ( left ), 𝑄 𝜙 ( middle ) and 𝑈 𝜙 ( right ) 𝐻 band imagesobtained with VLT/NACO using polarimetric differential imaging. No sig-nificant polarised signal is detected. The expected location of the companionis shown with a white circle. brary. The only requirement is to provide a snippet function to best_fit_tmp in order to read the template files. Either one or twofree parameters can be considered to find the best match: a flux scal-ing factor and, optionally, optical extinction. Two options are avail-able for the search of these optimal values: either a grid search (witha user-provided range) or a Nelder-Mead simplex algorithm, whichis faster and can also be constrained to an allowed range of values.The routine then returns a user-defined number of best-fit templatesminimizing the goodness-of-fit. As for mcmc_spec_sampling , thegoodness of fit takes into account spectral covariance and weights(i.e. it is given by Equation 2). Depending on the spectral resolu-tion of the template, it is either interpolated or convolved with thefilter(s) used for the observed spectrum. Figure 1c shows the final NACO image obtained for the 2017 dataset.Given the excellent and stable observing conditions, median-ADIwas sufficient to reveal a point source at a significant level at aseparation of 0 . (cid:48)(cid:48) ∼ 𝜎 detection.The 2017 detection motivated us to follow up the source in2018. The conditions during the 2018 visitor observing run weremediocre and the seeing (at 𝜆 = . (cid:48)(cid:48) . (cid:48)(cid:48) ∼ 𝑛 pc = 𝐻 band. We measured an SNR of 91 and 58 in the IRDIS 𝐾 𝐾 𝐾 𝐾
2. The companion was also recoveredby applying sPCA on the short non-saturated set of images, which Δ RA [mas]−40−200204060 Δ D E C [ m a s ] Figure 3.
Multi-epoch astrometry of the companion extracted using MCMC-NEGFC on the NACO 2017, NACO 2018 and IRDIS 2019 datasets, com-pared to predictions for a fixed background object based on the propermotion of CrA-9 measured by Gaia (Gaia Collaboration et al. 2018). A fixedbackground object can be rejected at a 5 𝜎 confidence level. consists of two individual cubes acquired at the beginning and endof the observation respectively (Table 2).In addition to the point source, the PCA-SDI image reveals atentative spiral pattern (Figure 1e). A possible primary arm extend-ing to the south of the disc appears to point towards the compan-ion candidate, as expected from hydro-dynamical simulations (e.g.Dong et al. 2015). We measure a signal-to-noise ratio (SNR) of ∼ 𝑛 pc = 𝑛 pc used. It is also tentativelyseen using PCA-ASDI, although self-subtraction of azimuthally ex-tended structures due to angular differential imaging may accountfor the differences between the two images. PCA-SDI and PCA-ASDI images obtained with different 𝑛 pc subtracted are providedin Appendix E. We also computed a standardized trajectory inten-sity mean map (STIM map; Pairet et al. 2019), using the residualcube after subtraction of the PCA-SDI model (as in Christiaenset al. 2019a). The STIM map does not reveal the spiral feature con-spicuously, as the inverse STIM map (i.e. obtained with oppositederotation angles) reveals pixels of similar intensity as the spiralfeature seen in the (regular) STIM map. Finally, we also applied themayonnaise algorithm on the core-saturated IRDIS dataset (Pairetet al. 2020). While the mayonnaise image may suggest some ex-tended disc signal to be present, the spiral pattern seen with IFS isnot recovered (Figure E2). New observations are thus required toconfirm the spiral pattern.Figure 2 shows the final 𝑄 𝜙 and 𝑈 𝜙 images. The central 𝑟 < . (cid:48)(cid:48) MNRAS , 1–22 (2021)
V. Christiaens et al.
Table 3.
Parameters inferred for the companion in our different datasets.Date Instrument Filter r PA contrast App. Mag a Abs. Mag b [mas] [ ◦ ] ( × − ) [mag] [mag]2017-08-26 NACO 𝐿 (cid:48) . ± . . ± . . ± . . ± .
18 10 . ± . 𝐿 (cid:48) . ± . . ± . . ± . . ± .
10 10 . ± . c — 699 . ± . . ± . . ± . d 𝐽 — — 8 . ± . . ± .
02 11 . ± . 𝐾 . ± . . ± . . ± . . ± .
03 10 . ± . 𝐾 . ± . . ± . . ± . . ± .
05 10 . ± . a Apparent magnitude. Uncertainties include the uncertainties on the contrast, the stellar flux photon noise and either theuncertainties on the stellar spectrum model (for the SPHERE measurements) or the uncertainties on the calibrated 𝐿 (cid:48) flux (for the NACO measurements). b Absolute magnitude after dereddening assuming 𝐴 𝑉 = . ± . c Reported uncertainties include both systematic and statistical (i.e. dispersion over all spectral channels) uncertainties. d Reported magnitude is integrated over the 2MASS 𝐽 band filter transmission curve. Only the 𝐽 -band filter is usedsince the IFS channels only partially cover the 𝑌 - and 𝐻 -band filter transmission curves. and combines images from a larger spectral bandwidth ( 𝑌 to 𝐻 bands instead of only 𝐻 ). We applied MCMC-NEGFC individually to each spectral channelof the IFS, the 𝐾 𝐾 𝐿 (cid:48) datasets. We show in Figure D1 three example corner plots obtainedby MCMC-NEGFC among our 43 ADI sequences, for the NACO2017, IFS 2019 (first spectral channel) and IRDIS 2019 ( 𝐾 Figure 3 shows the astrometric points retrieved by MCMC-NEGFCfor the companion candidate in the NACO 2017, NACO 2018 andSPHERE 2019 datasets. For the SPHERE 2019 epoch, we only con-sidered the IRDIS measurement given both the higher astrometricaccuracy and higher SNR of the point source than in the IFS im-ages. To make sure our astrometric extraction was accurate, we alsogot our dataset reduced by the SPHERE data center, and inferredconsistent astrometric estimates within 10% of our reported uncer-tainties. The IFS data were plagued by significant stripes hence notfurther considered in this work).We considered four sources of uncertainty: (i) the residualspeckle or background noise uncertainty captured by the varianceof the MCMC-NEGFC posterior distribution on 𝑟 and PA; (ii) astellar centering uncertainty conservatively assumed to be 0.1 pixelin the NACO and IRDIS datasets and 0.5 pixel in the IFS dataset(where the Strehl ratio was significantly lower); (iii) the systematicuncertainty on the plate scale, affecting the radial separation esti-mate; and (iv) the uncertainty on the PA of true north, affecting thePA estimate. We combined the different sources of uncertainty inquadrature for each parameter. For NACO, Launhardt et al. (2020)reported plate scale and true north measurements of 27 . ± . . ± . ◦ based on all their astrometric measure-ments between December 2015 and March 2018, respectively. Weadopted these values, but conservatively adopted an uncertainty of 0.5 ◦ for the PA of TN, to account for any difference between forthe different epochs of observations. We expect this uncertainty tobe conservative considering the consistent independent PA of TNestimate presented in Milli et al. (2017) based on 2016/09 data(0 . ◦ ± . ◦ ), and the maximum variation of 0.3 ◦ for the PA ofTN of NACO for all astrometric measurements within 2 years timereported in Chauvin et al. (2012). For IRDIS, we adopted the sys-tematic uncertainties quoted in Maire et al. (2016): a true north of − . ± .
08, and plate scales of 12.267 and 12.263 mas/px withthe 𝐾 𝐾 𝜎 confidence level that the object is a background star withnull proper motion. Instead, the measurements are consistent withnegligible orbital motion over the course of ∼ ∼
108 au. Our astrometric measurements areprovided in Table 3.
Figure 4 shows the contrast spectrum of the point source, i.e. the fluxratio with respect to the star at each wavelength. Our third modifica-tion to MCMC-NEGFC allows temporal variations throughout thedifferent observed sequences to be accounted for, hence enabling usto reach high precision on the estimated contrast of the point source.The accuracy of the contrast spectrum only depends on the residualspeckle and background noise level at the separation of the pointsource, which is captured by the variance of the posterior distribu-tion for the contrast (Figure D1). Although the SNR of the pointsource is higher in the IRDIS images obtained with the stellar-coresaturated dataset, it was not used to infer the contrast of the pointsource given the ignorance on the temporal variation of the stellarflux. Instead, the contrast in the 𝐾 𝐾 𝐾 𝐾 MNRAS000
Figure 4 shows the contrast spectrum of the point source, i.e. the fluxratio with respect to the star at each wavelength. Our third modifica-tion to MCMC-NEGFC allows temporal variations throughout thedifferent observed sequences to be accounted for, hence enabling usto reach high precision on the estimated contrast of the point source.The accuracy of the contrast spectrum only depends on the residualspeckle and background noise level at the separation of the pointsource, which is captured by the variance of the posterior distribu-tion for the contrast (Figure D1). Although the SNR of the pointsource is higher in the IRDIS images obtained with the stellar-coresaturated dataset, it was not used to infer the contrast of the pointsource given the ignorance on the temporal variation of the stellarflux. Instead, the contrast in the 𝐾 𝐾 𝐾 𝐾 MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? μ m)0.00080.00100.00120.0014 C on t r a s t IFSIFS - JIRDIS K1IRDIS K2NACO 2017NACO 2018
Figure 4.
Contrast spectrum of the companion extracted using MCMC-NEGFC on all unsaturated datasets. The overall positive slope with wavelength suggestsa redder companion than the central star. The two discrepant 𝐿 (cid:48) contrast estimates appear consistent with the variation in absolute stellar flux measured forCrA-9 in 2017 and 2018 (Section 3.1.2). model independent conclusions for the point source. Furthermore, itmay be used in future studies to re-estimate the companion’s (flux)spectrum if a higher-resolution calibrated stellar spectrum in theIR becomes available. Two features can be noted from the contrastspectrum:(i) The two NACO 𝐿 (cid:48) points are discrepant. To test a possiblebias related to the poorer quality of the NACO 2018 images, weran MCMC-NEGFC again with a variation of our third modifica-tion: instead of injecting the median unsaturated stellar PSF withvarying fluxes in the individual images of the cube, the injectionwas directly made with the corresponding stellar PSF, scaled tothe tested contrast. This led to a consistent contrast estimate for the2018 point, i.e. to the same level of discrepancy with the 2017 point.This suggests that the discrepancy may rather be due to variabilityof the primary star and/or the companion. The ∼
40% relative dif-ference between the 2017 and 2018 contrasts is consistent with thedifference in absolute flux obtained for CrA-9 using STD stars atthese 2 epochs: 23% larger and 19% smaller than the expected fluxat 𝐿 (cid:48) band based on the WISE W1 measurement of CrA-9 (Wrightet al. 2010). This suggests that the contrast discrepancy may beassigned entirely to the variability of the primary star, as it leads toa consistent 𝐿 (cid:48) flux for the point source in 2017 and 2018.(ii) We used the 2MASS 𝐽 -band filter transmission curve to inferthe contrast that would have been measured in that broad band filter(light blue point in Figure 4). Comparison of the 𝐽 -band contrast tothe 𝐾 𝐾 In order to infer the spectrum of the point source in contrast of thatof the star, an absolute calibration of the stellar flux measurementsis required first. This can be obtained through flux calibrators ob-served in similar conditions as the star of interest. Alternatively a reliable model spectrum for the star can be used. Given the im-possibility to obtain a good absolute calibration for the SPHEREstellar flux measurements (Section 3.1.3), we opted for the secondoption. Although past studies have estimated the spectral type andeffective temperature of CrA-9, the spectrum of the star is currentlypoorly sampled at IR wavelengths. Therefore, we used specfit incombination with the BT-SETTL grid of models to infer the mostlikely SED for CrA-9 in the 0.9–4.0 𝜇 m range. BT-SETTL modelsconsider a parameter-free cloud prescription to account for dust for-mation, coagulation and settling (Allard et al. 2012; Allard 2014).We only considered measurements made by (1) Gaia (DR2) in its 𝐺 𝑅𝑃 filter (Gaia Collaboration et al. 2016, 2018), (2) 2MASS inthe JHK bands (Cutri et al. 2003), and (3) WISE at 3.3 and 4.6 𝜇 m(W1 and W2 bands; Wright et al. 2010). These instruments havethe smallest reported photometric uncertainties in their respectivewavelength range. We did not extend the wavelength range to avoidthe need to add extra components to the model to account for eitheraccretion luminosity or disc emission, and to be less affected by thepoorer knowledge of the extinction law at short wavelength. Includ-ing more photometric points at either shorter or longer wavelengthswould involve more biases and likely lead to a poorer model in therange of interest.T-Tauri stars are known to show significant variability overtime due to chromospheric activity and/or accretion (e.g. Bouvieret al. 2004; Hartmann et al. 2016, and references therein). Our 𝐿 (cid:48) stellar flux measurements calibrated using standard stars suggestthat CrA-9 is no exception (Section 3.1.2). This can also be seenfrom the vertical scatter of points at different optical to near-IRwavelengths in the left panel of Figure 5. In order to take into accountthe effect of variability, we considered the reported photometricuncertainties of the 𝐽𝐻𝐾 𝑇 e ≈ ±
150 K, log ( 𝑔 ) ≈ . ± . 𝐴 𝑉 ≈ . ± . 𝑅 𝑉 in our model,to account for possibly different grain sizes in the line of sight than MNRAS , 1–22 (2021) V. Christiaens et al. Wavelength ( μ m)10 −18 −17 −16 −15 −14 −13 λ F λ ( W m − ) SED of CrA-9Points used for the fitBT-SETTL model (T e =3630K; A V =1.7mag) 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5Wavelength ( μ m)0.00.51.01.52.02.53.03.5 1e−13 SED of CrA-9Favoured BT-SETTL model (sampled as observations)BT-SETTL models drawn from posterior distribution Figure 5. ( left ) Spectral energy distribution of CrA-9 ( black points) and favoured BT-SETTL model as inferred by specfit ( dark green curve). The lightgreen points are the ones used for the fit. ( right ) Zoom on the wavelength range relevant for our observations. The dark green points indicate the flux of thefavoured BT-SETTL model at the sampling of our IFS, IRDIS and NACO observations. The blue curves are randomly sampled models from the inferredposterior distribution (Figure F1), which are used to estimate the uncertainties on the model spectrum of the star at each wavelength. in the diffuse ISM (e.g. Weingartner & Draine 2001; Calvet et al.2004). Figure 5 shows the model favoured by specfit in green.The right panel of Figure 5 provides a zoom on the wavelengthrange of interest and shows 60 sample models from the posteriordistribution (in light blue). The associated corner plot is shown inFigure F1. For each parameter, the quoted uncertainties correspondto the 34th and 66th percentile of the posterior distribution. Wefind an effective temperature 𝑇 e ≈ + − K, a surface gravitylog ( 𝑔 ) ≈ . ± .
2, a photometric radius 𝑅 ≈ . + . − . 𝑅 (cid:12) , anoptical extinction 𝐴 𝑉 ≈ . + . − . mag and a broad constraint on theratio of total-to-selective optical extinction that is consistent withthe diffuse ISM value.We drew 1000 sample spectra from our posterior distribu-tion in order to estimate uncertainties on the stellar spectrum atthe wavelengths of our observations. We fitted a Gaussian profileto the distribution of sample values at each relevant wavelengths,and considered the standard deviation of each Gaussian fit as theuncertainty in stellar flux, for error propagation to the companioncandidate spectrum.We note that the 𝐿 (cid:48) flux of the posterior BT-SETTL models(1 . ± . × − W m − 𝜇 m − ) is consistent with the aver-age of the two absolute 𝐿 (cid:48) fluxes estimated in the 2017 and 2018datasets through standard stars (1 . ± . × − W m − 𝜇 m − ;Section 3.1.2). The spectrum of the point source is obtained by multiplying our con-trast spectrum (Figure 4) to (i) the favoured BT-SETTL model of thestar (Figure 5) for the SPHERE points, after the model is convolvedwith the relevant filters used in the SPHERE observation; and (ii)the absolute 𝐿 (cid:48) flux measurements for the star for the 2017 and 2018NACO points. For the rest of our analysis, we adopt the weightedaverage of the two 𝐿 (cid:48) flux values, which are consistent with eachother. The spectrum of the faint point source is shown with blackpoints in Figure 6, with the two individual 𝐿 (cid:48) measurements for2017 and 2018 shown with gray points . Our final uncertainties onthe spectrum include three contributions combined in quadrature: (1) the uncertainties on the contrast inferred for the point source byNEGFC-MCMC (Section 4.2.2); (2) the photon noise on the mea-sured flux of the star in each dataset; and (3) the uncertainty on thebest-fit BT-SETTL model of the star (Section 4.3). The latter domi-nate the uncertainty budget, being an order of magnitude larger thanthe uncertainties on the contrast and the stellar photon noise at allwavelengths. We report in Table 3 the apparent magnitudes inferredfor the point source in the 𝐽 , 𝐾 𝐾 𝐿 (cid:48) filters considering theabove uncertainties. The most significant feature in the spectrum of the point source isthe broad H O absorption band spanning ∼ 𝜇 m (Auman1967), which has been observed in a number of young low-masscompanions (e.g. Bonnefoy et al. 2014). We compute the H Oand Na indices introduced in Allers et al. (2007) to estimate thespectral type and gravity of M5 to L5 dwarfs, respectively. TheH O index is insensitive to gravity and the Na index is relativelyinsensitive to spectral type, contrary to most other spectral indiceswhich are sensitive to both (Allers et al. 2007; Bonnefoy et al. 2014).We measured an H O index of (cid:104) 𝐹 𝜆 = . − . (cid:105)/(cid:104) 𝐹 𝜆 = . − . (cid:105) = . ± .
007 and a Na index of (cid:104) 𝐹 𝜆 = . − . (cid:105)/(cid:104) 𝐹 𝜆 = . − . (cid:105) = . ± . ± ∼ 𝐴 𝑉 = ∼ 𝐿 (cid:48) measurements, we only considered theIFS+IRDIS spectrum of CrA-9 and compared it to all MSL tem- MNRAS000
007 and a Na index of (cid:104) 𝐹 𝜆 = . − . (cid:105)/(cid:104) 𝐹 𝜆 = . − . (cid:105) = . ± . ± ∼ 𝐴 𝑉 = ∼ 𝐿 (cid:48) measurements, we only considered theIFS+IRDIS spectrum of CrA-9 and compared it to all MSL tem- MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? μ m)0.00.51.01.52.02.5 λ F λ ( W m − ) γ ? Br γ BT-SETTL: T e =3144K, R b =0.57 R J ( Δ AIC = 0)BT-DUSTY: T e =3074K, R b =0.60 R J ( Δ AIC = -7)BT-DUSTY: T e =3059K, R b =0.60 R J (with Br γ ; Δ AIC = -15)BT-DUSTY: T e =3109K, R b =0.59 R J (IFS+IRDIS only)DRIFT-PHOENIX: T e =3000K, R b =0.60 R J ( Δ AIC = 62)
Figure 6.
Measured spectrum of CrA-9B/b ( black points ) compared to the favoured BT-SETTL ( blue ), BT-DUSTY ( green ) and DRIFT-PHOENIX ( yellow )models returned by specfit . The results are obtained with 𝑅 𝑉 set to 3.1, although the curves obtained when 𝑅 𝑉 is set as a free parameter are visuallyidentical. The dashed green curve considers the luminosity of a Br 𝛾 emission line as free parameter. The dotted green curve is obtained considering only theIFS+IRDIS points for the fit. The measurements considered for the fit are shown in black . At 3.8 𝜇 m, we consider the weighted average of the NACO 2017 and2018 individual measurements shown in gray . Including a Br 𝛾 emission line reproduces better all measurements, including the 𝐾 plates (both field and young) with a wavelength range covering0.9–2.3 𝜇 m (i.e. 326 templates in total). We let both scaling andextinction as free parameters in our fit, and used the simplex searchmode of the best_fit_tmp routine. Figure 7 shows the best threefits ( 𝜒 𝑟 ∼ . . 𝐴 𝑉 ∼ 𝛽 Pic moving group (Gagné et al.2015). These results are consistent with the spectral type and low-gravity inferred using the H O index and Na index.
To retrieve physical parameters for the companion candidate, wesubsequently ran specfit with the BT-SETTL, BT-DUSTY andDRIFT-PHOENIX grids of models. We considered four free param-eters for the BT-SETTL and BT-DUSTY models ( 𝑇 e , log ( 𝑔 ) , 𝑅 and 𝐴 𝑉 ), and five free parameters, adding the metallicity log ( 𝑍 / 𝑍 (cid:12) ) , for the DRIFT-PHOENIX models. All three grids are based onthe PHOENIX atmosphere models (Hauschildt 1992). Comparedto the BT-SETTL grid, BT-DUSTY models consider the maxi-mum amount of dust allowed in equilibrium with the gas phase,without a cloud prescription (Allard et al. 2011, 2012). The DRIFT-PHOENIX models use a kinetic approach to model dust forma-tion, growth, settling and advection (Woitke & Helling 2003, 2004;Helling & Woitke 2006; Helling et al. 2008).For all grids, we considered uniform priors on all parameters. Inthe case of the BT-SETTL grid, we considered all available modelswith 𝑇 𝑒 ∈ [ , ] K in steps of 100 K, and log ( 𝑔 ) ∈ [ . , . ] in steps of 0.5dex. For the BT-DUSTY grid, we considered allmodels with 𝑇 𝑒 ∈ [ , ] K and log ( 𝑔 ) ∈ [ . , . ] , owingto incomplete wavelength coverage for models at low temperatures.Similar results were obtained with either grids. We considered allavailable DRIFT-PHOENIX models with 𝑇 𝑒 ∈ [ , ] K insteps of 100 K, log ( 𝑔 ) ∈ [ . , . ] in steps of 0.5dex, and metallicitylog ( 𝑍 / 𝑍 (cid:12) ) ∈ [− . , , . ] . MNRAS , 1–22 (2021) V. Christiaens et al. μ m)0.500.751.001.251.501.752.002.252.50 λ F λ ( W m − ) A V =1.7mag)Template A V =2.1mag)Template A V =2.4mag)Measured spectrum of CrA-9B/b Figure 7.
SPHERE (IFS+IRDIS) spectrum of CrA-9 B/b ( black points) compared to the three best-fit templates of the Montreal Spectral Library: 2MASSJ05071137+1430013 B ( 𝜒 𝑟 ∼ .
7; Gagné et al. 2015), 2MASS J03390160-2434059 ( 𝜒 𝑟 ∼ .
7; Gagné et al. 2015) and GU Psc A ( 𝜒 𝑟 ∼ .
9; Naud et al.2014). The fit considered two free parameters: a scaling factor and optical extinction (whose best-fit value is provided in parenthesis).
The favoured models for the different grids are shown in Fig-ure 6, and the favoured parameters are provided in Table 4. For allgrids of models, the inferred effective temperature (3000–3200K)is consistent with a mid-M dwarf. However, the photometric ra-dius required to account for the faint measured flux is only ∼ 𝜎 ), except forthe 𝐾 ∼ 𝜎 outlier andthe 1.09 𝜇 m IFS channel (2–2.5 𝜎 ). The best-fit DRIFT-PHOENIXmodels hit the upper bound of the grid (3000 K), which accountsfor the poorer visual fit (Figure 6) and underestimated uncertain-ties on the different parameters (Table 4). The favoured value ofextinction ( 𝐴 𝑉 ≈ . ± . 𝐴 𝑉 ≈ . + . − . mag; Section 4.3). The surface gravity valuesappear slightly larger for the BT-SETTL and BT-DUSTY modelsthan expected from a 1–2 Myr-old mid-M dwarf ( ∼ 𝛾 emission line affecting the 𝐾 Δ AIC between its AIC and the AIC ob-tained with the grid of BT-SETTL models using 4 free parameters( 𝑇 𝑒 , log ( 𝑔 ) , 𝑅 and 𝐴 𝑉 ), and reported that difference in Table 4. Thelower the value of AIC (hence Δ AIC), the higher the likelihood.Our first test consisted in fixing the photometric radius to 1.8 𝑅 Jup , letting the other parameters free. The best-fit BT-SETTL andDRIFT-PHOENIX models obtained with such constraint are shownin Figure G2. These models are poor fits to the data, with Δ AIC val-ues larger than 2000. This suggests that the solution found withoutconstraint ( 𝑇 e ∼ 𝑅 𝑝 ∼ . 𝑅 Jup ) results indeedfrom the lack of good solutions with larger photometric radii.Next we tested the addition of another free parameter to modeldust extinction: 𝑅 𝑉 . Since the system is young, a bound companionmay be surrounded by its own disc of gas and dust. Dust growth indiscs may affect the extinction law, possibly leading to a differentvalue of 𝑅 𝑉 than in the diffuse ISM (e.g. Weingartner & Draine2001; Calvet et al. 2004). We compare in Figure 8 the posteriordistributions of the BT-DUSTY model parameters yielded by spec-fit when 𝑅 𝑉 is set as a free parameter or not (similar results areobtained with BT-SETTL). The favoured values of 𝑇 𝑒 , log ( 𝑔 ) and 𝑅 are consistent with the results obtained when 𝑅 𝑉 is fixed (to 3.1).The only difference is a significant degeneracy between the valuesof 𝐴 𝑉 and 𝑅 𝑉 , with a range of values for ( 𝐴 𝑉 , 𝑅 𝑉 ) pairs all leadingto similar quality fits. Considering or not 𝑅 𝑉 as a free parameter ledto similar AIC values for either BT-SETTL or BT-DUSTY models(Table 4), and none of the favoured models is able to reproduce the 𝐾 Δ AIC <
10 betweenall models which fixed 𝑅 𝑉 or not, there is no significant supportfor one of these models over the others (e.g. Burnham & Anderson2002). MNRAS000
10 betweenall models which fixed 𝑅 𝑉 or not, there is no significant supportfor one of these models over the others (e.g. Burnham & Anderson2002). MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? Table 4.
Physical properties of CrA-9B/b inferred from models with high support retrieved by specfit for the different atmospheric grids.Parameter From spectral BT-SETTL BT-DUSTY DRIFT-PHOENIX BT-DUSTY BT-DUSTY BT-DUSTYindices (fixed 𝑅 𝑉 ) (fixed 𝑅 𝑉 ) (fixed 𝑅 𝑉 ) (free 𝑅 𝑉 ) (free 𝐿 Br 𝛾 ) (SPHERE only)Spectral type M5.5 ± ( a ) — — — — — — 𝑇 e [K] 2910 ± ( b ) + − ±
43 3000 + − ( e ) + − + − + − Log( 𝑔 ) 3 . ± . ( c ) . + . − . . + . − . . ± . . + . − . . + . − . . + . − . Radius [ 𝑅 𝐽 ] 7 . ± . ( d ) . ± .
01 0 . + . − . . ± .
01 0 . + . − . . ± .
01 0 . ± . 𝐴 V [mag] — 2 . ± . . ± . . + . − . . + . − . . ± . . ± . 𝑅 V — (3.1) (3.1) (3.1) 1 . + . − . (3.1) 1 . + . − . log ( 𝑍 / 𝑍 (cid:12) ) — — — − . + . − . ( e ) — — —log ( 𝐿 Br 𝛾 𝐿 (cid:12) ) — — — — — − . + . − . — Δ AIC — 0 − − −
15 — a Based on the H O index and empirical relation derived in Allers et al. (2007). b Based on the empirical relation to convert from spectral type to effective temperature in Herczeg & Hillenbrand (2014). c Considering the Na index value, an age of 2 Myr (since the Na index is consistent with that measured for Cha I members; Allers et al. 2007),and the isochrones of either Tognelli et al. (2011) or Baraffe et al. (2015) to estimate a log ( 𝑔 ) value. d Considering 𝑇 e = ± e Parameter hits bound of the grid; reported uncertainties for all parameters of this model are to be considered lower limits.
Given the possibility for the companion to be variable likethe host star, it is unclear whether its 𝐿 (cid:48) flux is the same at theepoch of the SPHERE dataset and at the epochs of the NACOdatasets. Furthermore, our absolute flux calibration of the SPHEREmeasurements is based on CrA-9 which is known to be variable,and could hence lead to a shift with respect to the 𝐿 (cid:48) measurement.Therefore we also ran specfit on the IFS+IRDIS points only. Inthat case, we found that the favoured BT-SETTL and BT-DUSTYmodels are consistent with the ones obtained when including the 𝐿 (cid:48) point. This suggests that the puzzling parameter values found aboveand the impossibility for the models to reproduce the 𝐾 𝛾 emission line ( 𝜆 = . 𝜇 m) as a free parameter, in an at-tempt to account for the observed discrepancy between the 𝐾 𝛾 emission, butfixing 𝑅 𝑉 to 3.1 and including the weighted-average NACO point.specfit assumes a Gaussian profile for injected lines. To limit theline injection process to a single free parameter – its flux, we setthe full width at 10% height to 100 km s − , based on the observedH recombination line width of the PDS 70 and Delorme 1 (AB) bplanets (Haffert et al. 2019; Eriksson et al. 2020). Given the lowresolution of our spectrum and the 𝐾 𝛾 emission line enable us to reproduce the observed 𝐾 𝛾 luminosity log ( 𝐿 Br 𝛾 𝐿 (cid:12) ) ≈ − . ± .
1, and lead tosimilar physical parameters as estimated without including that freeparameter (Figure 6; Table 4). The values of Δ AIC for these types ofmodel ( − −
13 for BT-SETTL and BT-DUSTY, respectively)suggests marginal support in favour of including the Br 𝛾 luminosityas a free parameter. The favoured BT-DUSTY model including Br 𝛾 line emission is shown with a green dashed line in Figure 6, and thecorresponding corner plot is shown in Figure 8. In this section, we discuss several possible interpretations for ourresults, their respective likelihood, and the kind of observationsrequired to discrimate our two leading hypotheses.
A background star scenario would explain why the BT-SETTL andBT-DUSTY models favoured by specfit have an anomalously smallphotometric radius ( ∼ . 𝑅 𝐽 at the distance of CrA-9), and a largersurface gravity (log ( 𝑔 ) ∼ .
5) than expected for an object of theage of CrA-9 (log ( 𝑔 ) ∼ . ∼ 𝑅 𝐽 (Baraffeet al. 2015), hence lying at a distance several times that of CrA-9.However, the three-epoch astrometric positions we measured areconsistent with a companion co-moving with CrA-9 (Figure 3).We can consequently reject a background object with negligibleproper motion with a 5 𝜎 confidence. Furthermore, it is unlikelythat a background star several times further than CrA-9 had a highproper motion perfectly matching the significant foreground mo-tion of CrA-9. A high-proper motion background star scenario isalso inconsistent with the measured 1.13 𝜇 m-Na spectral index (Sec-tion 4.4.2). The value of 1 . ± .
007 is consistent with a youngvery low-gravity object and rejects at a 5 𝜎 confidence the samegravity as field dwarfs.The Galactic latitude of CrA-9 ( 𝑏 ≈ − ◦ ) further makes abackground star several times the distance to CrA-9 unlikely. Ourestimated probability of the point source being a background staris (cid:46) ∼ . (cid:48)(cid:48) 𝐿 (cid:48) band. We estimated this by considering a spatialhomogeneous Poisson point process with the rate set by the numberdensity of objects brighter than 𝐿 (cid:48) = . MNRAS , 1–22 (2021) V. Christiaens et al. T e = +43−43 K . . . . l og ( J ) log( J ) = +0.2−0.4 . . . . R b R b = +0.01−0.02 R J T e . . . . A V . . . . log( J ) .
54 0 .
57 0 .
60 0 . R b .
50 1 .
75 2 .
00 2 . A V A V = +0.1−0.1 mag T e = +56−36 K . . . . l og ( J ) log( J ) = +0.2−0.4 . . . . R b R b = +0.01−0.02 R J . . . . A V A V = +0.4−0.2 mag T e R V . . . . log( J ) .
54 0 .
57 0 .
60 0 . R b . . . . A V R V R V = +3.0−0.2 T e = +49−36 K . . . . l og ( J ) log( J ) = +0.2−0.4 . . . . . R b R b = +0.01−0.01 R J . . . . L B r γ ( L ⊙ ) log( L BrG L ⊙ ) = −5.89 +0.06−0.10 T e . . . . . A V . . . . log( J ) . . . . . R b . . . . L Brγ ( L ⊙ ) . . . . . A V A V = +0.1−0.1 mag Figure 8.
Corner plots retrieved by specfit using the BT-DUSTY modelsto reproduce the spectrum of CrA-9 B/b, when ( top ) fixing 𝑅 𝑉 to thediffuse ISM value, ( middle ) setting 𝑅 𝑉 as a free parameter, and ( bottom )considering the Br 𝛾 luminosity as a free parameter. Considering the above, the evidence is in favour of a boundcompanion. We thereafter refer to the point source as companion . Considering a distance of 𝑑 = . ± . . ± . 𝑅 𝑉 = . . ± .
10 mag in the 𝐽 band, 10 . ± .
07 mag with the 𝐾 𝐾 . ± .
19 mag or 10 . ± .
11 mag in the 𝐿 (cid:48) band, based on the 2017 and 2018 epochs data respectively (Table 3).Considering an age of 1–2 Myr, these absolute magnitudes wouldcorrespond to masses of ∼ 𝑀 𝐽 according to either the CONDevolutionary models (Baraffe et al. 2003) or the hot-start modelspresented in Spiegel & Burrows (2012). However the 𝐽 − 𝐿 (cid:48) colorof 1.3–1.7 mag appears bluer than predicted by those models ( > . 𝑅 𝐽 is also significantlysmaller than expected for a gas-dominated planet, and the estimatedeffective temperature of the companion (3000–3200 K) appears todefy even the most optimistic hot-start models for a 10 𝑀 𝐽 planet(e.g. Baraffe et al. 2003; Spiegel & Burrows 2012; Mordasini et al.2012). The companion also shows a significantly bluer spectrumthan PDS 70 b (Müller et al. 2018; Christiaens et al. 2019b), althoughthis might be due to PDS 70 b being more significantly enshroudedby dust (Christiaens et al. 2019b; Wang et al. 2020; Stolker et al.2020b).The arguments above a priori suggest that CrA-9B/b is un-likely to be a protoplanet. Nonetheless, an alternative explanationfor the observed spectrum is that it traces the fraction of a proto-planet photosphere that is heated by an accretion shock (e.g. Zhu2015; Aoyama et al. 2020). This scenario would be compatible witha significant Br 𝛾 emission line (2.166 𝜇 m), possibly required to ac-count for the observed 𝐾 𝛾 luminosity log ( 𝐿 Br 𝛾 / 𝐿 (cid:12) ) ≈ − . ± . (cid:38) 𝜎 difference between the favouredBT-DUSTY and BT-SETTL models (without Br 𝛾 ) and the mea-sured 𝐾 𝛾 line is a known tracer of gasaccretion for classical T-Tauri stars (Muzerolle et al. 1998; Calvetet al. 2004), and also possibly for protoplanets (Aoyama et al. 2020).According to the protoplanet accretion shock models presented inAoyama et al. (2020), a log ( 𝐿 Br 𝛾 / 𝐿 (cid:12) ) ≈ − . ∼ − 𝑀 𝐽 yr − . This estimate is signifi-cantly larger than inferred for the PDS 70 planets (Haffert et al.2019; Hashimoto et al. 2020), although the latter may be at a moreadvanced stage considering they have already cleared a wide anddeep gap (e.g. Hashimoto et al. 2011; Dong et al. 2012; Keppleret al. 2020). The inferred mass accretion rate is similar to whathas been predicted from magneto-hydrodynamical simulations ofaccreting sub-Jovian planets embedded in the protoplanetary disc(e.g. Gressel et al. 2013). This would require the companion tostill lie within the gaseous component of the protoplanetary disc.Although the estimated dust disc radius ( ∼ . (cid:48)(cid:48)
39; Cazzoletti et al.2019) is smaller than the separation of the companion, there is nocurrent constraint on the extent of the gas disc which is likely toextend to larger separations than the mm-size dust, hence possiblyup to the separation of the companion.If we are indeed witnessing an accreting protoplanet, a sig-nificant fraction of the observed luminosity would be expected tocome from accretion luminosity (potentially larger than the photo-spheric contribution; e.g. Mordasini et al. 2017; Marleau et al. 2017,2019b). Whether the accretion stream shocks on the circumplan-
MNRAS000
MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? etary disc or directly onto the photosphere of the protoplanet, thesurface is expected to reach (cid:38) (cid:38) ∼ − 𝑀 𝐽 yr − (based on 𝐿 Br 𝛾 and the models in Aoyama et al.2020), the observed blue slope in the 1.0–3.8 𝜇 m range may be com-patible with a subset of the models presented in Zhu (2015) withdifferent combinations of inner truncation radius for the circum-planetary disc and filling factors between 1–10% (i.e. the fractionof the surface covered by accretion shocks). Finally, the above sce-nario would account for the small photometric radius inferred in thiswork; a 0.6 𝑅 𝐽 photometric radius makes for 2–10% of the area ofa 1.6–5 𝑅 𝐽 protoplanet predicted from cold- to warm-start modelsin Aoyama et al. (2020).A potential caveat of the above hypothesis is that one may alsoexpect to observe Pa 𝛽 (1.282 𝜇 m) and Pa 𝛾 (1.094 𝜇 m) recombina-tion lines in emission (e.g. Aoyama et al. 2020). A 2.5 𝜎 excess canbe seen at ∼ . 𝜇 m but no excess is seen at ∼ 𝜇 m. However,our conclusions on the H recombination lines of the companion arelimited by our ignorance of the emission lines affecting the star.Indeed, the spectrum of the companion is obtained in contrast of amodel SED for the star which does not include its emission lines (i.e.only the physics included in the BT-SETTL/BT-DUSTY models).The star is known to have a significant H 𝛼 emission line though(Romero et al. 2012), which suggests it also harbours other H re-combination lines such as Pa 𝛽 , Pa 𝛾 or Br 𝛾 . Not including those inthe model SED of the star leads to lower limits on the inferred flux ofthose lines in the companion spectrum. It is nevertheless possible todiscuss the relative line ratios between the star and the companion.If the excesses at ∼ 𝜇 m and in the 𝐾 𝛾 and Br 𝛾 emission lines, respectively, the lack of excess at ∼ 𝜇 m suggests that both the Pa 𝛾 /Pa 𝛽 and Br 𝛾 /Pa 𝛽 line ratiosare larger for the companion than for the star. This translates toPa 𝛾 𝑏 /Pa 𝛽 𝑏 (cid:38) . 𝛾 𝑏 /Pa 𝛽 𝑏 (cid:38) .
3, where we considered thereported luminosity and mass accretion rate of CrA-9 to estimate itsexpected line ratios using the models in Edwards et al. (2013). Theformer requirement may be met by either an accreting T-Tauri star orprotoplanet for reasonable assumptions on the shock velocity and Hnumber density (Edwards et al. 2013; Aoyama et al. 2020), but thelatter condition appears more difficult to reconcile with predictionsfor an accreting protoplanet (Br 𝛾 𝑏 /Pa 𝛽 𝑏 ≈ ∼
10 yearsearlier (H 𝛼 measurement presented in Romero et al. 2012) than our2019 SPHERE observations; and that the filling factors assumed forline ratio predictions were ∼
10% (Edwards et al. 2013) and ∼ The recent examples of FW Tau C (Kraus et al. 2014; Wu & Shee-han 2017) and CS Cha B (Ginski et al. 2018; Haffert et al. 2020)show that highly extinct stellar binary companions may mimic thesignal of a planetary-mass companion, in particular when the onlyavailable information is a near-IR flux and/or colour. In the case ofthe faint companion around CrA-9, the optical extinction derivedby the spectral analysis appears moderate ( 𝐴 𝑉 ∼ . 𝑀 (cid:12) M-dwarf for an age of 1–2 Myr (e.g. Baraffe et al.2015). Considering the 𝐿 Br 𝛾 – 𝐿 acc empirical relations calibrated onlow-mass T-Tauri stars (Muzerolle et al. 1998), the Br 𝛾 luminosityof log ( 𝐿 Br 𝛾 / 𝐿 (cid:12) ) ≈ − . ( 𝐿 acc / 𝐿 (cid:12) ) ≈ − .
0. This converts toa mass accretion rate of ∼ × − 𝑀 (cid:12) yr − for a young 0.15- 𝑀 (cid:12) M dwarf. This is similar to the average accretion rate of low-massT-Tauri stars (6 × − 𝑀 (cid:12) yr − ; Gullbring et al. 1998; Calvet et al.2004), and to the mass accretion rate expected for CrA-9 based onits measured H 𝛼 luminosity (2 . × − 𝑀 (cid:12) yr − ). A significantBr 𝛾 emission line has also been observed for HK Tau B (McCabeet al. 2011), which is another property shared with the companionseen around CrA-9.The obscured binary scenario may appear slightly at odds withthe lack of disc signal from that location. Our NACO PDI obser-vations did not detect any polarised signal at the location of thecompanion, despite good and stable observing conditions and a ∼
1h integration (Appendix 3.1.1). However, it is unclear whether theobservation was sensitive enough to detect the polarised fractionof the faint flux received from the companion. No sub-mm contin-uum emission was detected at the location of the companion either(Cazzoletti et al. 2019). Considering their RMS noise and a dusttemperature of 20 K, this corresponds to a 3 𝜎 upper limit on thedust mass of ∼ . 𝑀 ⊕ (or ∼ . 𝑀 𝐽 total disc mass assuming astandard 100:1 gas-to-dust ratio). The lack of resolved near-IR emis-sion around the point source translates into an upper limit of ∼
12 au(considering a 1.3FWHM resolution power achieved with IRDIS)for the extent of the circum-secondary disc, which corresponds to afraction of the Hill radius of a 0.1 𝑀 (cid:12) companion orbiting at (cid:38) (cid:38)
44 au Hill radius). This limit may be consistentwith the expected size of a small edge-on scattering disc producingall the flux observed from the companion. To obtain an order of mag-nitude estimate we consider the example of the well-characterisedHK Tau B edge-on circum-secondary disc. The disc of HK Tau B
MNRAS , 1–22 (2021) V. Christiaens et al. has a radius of 104 au which leads to a ∼
10 times underluminousM2 dwarf spectrum; we would thus expect a ∼
10 times smaller discradius for the ∼ ∼
10 au.A similar caveat affects the obscured accreting binary and ac-creting protoplanet scenarios: the lack of significant Pa 𝛽 in the IFSspectrum. The latter is observed to be brighter than the Br 𝛾 line formost CTTS (e.g. Calvet et al. 2004; Edwards et al. 2013). In thecase of an obscured binary, it is also surprising that such a signifi-cant Br 𝛾 emission line would not be as affected by obscuration asthe continuum flux. This may require an unlikely viewing geome-try. Alternatively, this leaves the door open to the possibility of themismatch between the 𝐾 𝛾 lineemission. Either the assumptions made in these models not applyingto the case of CrA-9 B/b (e.g. assumptions on atmosphere micro-physics, opacity sources, and/or metallicity), or another emissionline may be affecting the 𝐾 ( 𝑔 ) (cid:38) . ( 𝑔 ) ≈ .
7; Tognelliet al. 2011; Baraffe et al. 2015). A similar conclusion was reachedupon characterisation of the 𝐻 + 𝐾 spectrum of the young accret-ing M-dwarf companion HD 142527 B (log ( 𝑔 ) (cid:38) .
5; Christiaenset al. 2018). Could the large surface gravity values inferred for bothCrA-9 B/b and HD 142527 B be due to the BT-SETTL and BT-DUSTY models not including the effect of magnetic fields? Wheninferring physical parameters of low-mass T-Tauri stars from theiroptical/near-IR spectrum, not considering the effect of their strongmagnetic field (e.g. Johns-Krull et al. 1999; Johns-Krull 2007) isknown to lead to significant discrepancies in effective temperatureand surface gravity estimates (see e.g. Sokal et al. 2018, for the caseof TW Hya). However, to our knowledge this effect has not beeninvestigated on near-IR continuum shape, it is thus unclear whetherthis could account for the large inferred surface gravity.Finally, it is worth noting that if the point source is an obscuredM-dwarf, it could either trace a bound companion or a stellar fly-by. Either scenarios would be compatible with the tentative spiralsseen in the IFS image and the presence of a possibly large amountof obscuring material (e.g. Zhu et al. 2015; Cuello et al. 2020).Moreover, both scenarios are common outcomes in hydrodynamicalsimulations of star formation in relatively dense environments (e.g.Bate 2018).
The following types of follow-up observations would be the mostuseful to constrain the nature of the companion:(i) VLT/MUSE observations would set independent constraintson the accretion luminosity and mass accretion rate of the compan-ion and potentially provide a spectrum at visible wavelengths for thecompanion (e.g. Haffert et al. 2019; Hashimoto et al. 2020; Haffertet al. 2020).(ii) Either VLT/GRAVITY or Keck/KPIC observations couldprovide a medium resolution spectrum of the companion in the 𝐾 band, hence confirm whether a Br 𝛾 emission line is present(Gravity Collaboration et al. 2019; Mawet et al. 2016). Furthermore,GRAVITY could also constrain the size of the emitting region, as itwas recently done for PDS 70 b (Wang et al. 2021). In particular, it would be able to resolve an edge-on circum-secondary disc as smallas ∼ CO to probe the gaseouscomponent. The latter could provide an estimate of the mass of thecompanion based on the disc kinematics (e.g. Wu & Sheehan 2017,for the case of FW Tau C). These observations may also confirmsubstructures in the circumprimary disc (such as the tentative IRspirals), which could then be used to derive independent mass es-timates on the companion through hydro-dynamical modeling (e.g.Price et al. 2018; Calcino et al. 2020).(v) A combination of high-precision radial velocity and Gaia as-trometric measurements for CrA-9, together with our direct imagingconstraints, would enable to set independent constraints on the massof the companion (e.g. Brandt et al. 2019, 2020).Each of these observations could then be compared to the predic-tions from the scenarios presented in Sections 5.1.2 & 5.1.3, inorder to constrain the nature of the companion.
In this work, we developed the following methods:(i) We implemented new reduction pipelines for non-coronagraphic data obtained with VLT instruments NACO, IRDISand IFS. Each of these pipelines makes use of vip routines, whilethe latter two also make use of esorex recipes.(ii) We adapted the negfc module of the open-source packagevip (Gomez Gonzalez et al. 2017) in order to refine the extractionof the astrometry and contrast of directly imaged companions.(iii) We implemented specfit, a new module for the spectralcharacterisation of both stellar and substellar objects in a Bayesianframework, added it to vip, and used it to infer the most likely phys-ical parameters for both the star and the faint companion discoveredin this work.We have applied the methods listed above to analyse our dataon CrA-9. Our scientific results are summarised as follow:(i) We observed the T-Tauri star CrA-9, a known transition discwith a dust cavity, and detected a faint point source at 0 . (cid:48)(cid:48) 𝜎 confidence level, and isconsistent with a bound companion.(v) We determined that the companion was 7.1–7.9 mag fainterthan the star in the 1.0–3.8 𝜇 m wavelength range, leading to absolute MNRAS000
In this work, we developed the following methods:(i) We implemented new reduction pipelines for non-coronagraphic data obtained with VLT instruments NACO, IRDISand IFS. Each of these pipelines makes use of vip routines, whilethe latter two also make use of esorex recipes.(ii) We adapted the negfc module of the open-source packagevip (Gomez Gonzalez et al. 2017) in order to refine the extractionof the astrometry and contrast of directly imaged companions.(iii) We implemented specfit, a new module for the spectralcharacterisation of both stellar and substellar objects in a Bayesianframework, added it to vip, and used it to infer the most likely phys-ical parameters for both the star and the faint companion discoveredin this work.We have applied the methods listed above to analyse our dataon CrA-9. Our scientific results are summarised as follow:(i) We observed the T-Tauri star CrA-9, a known transition discwith a dust cavity, and detected a faint point source at 0 . (cid:48)(cid:48) 𝜎 confidence level, and isconsistent with a bound companion.(v) We determined that the companion was 7.1–7.9 mag fainterthan the star in the 1.0–3.8 𝜇 m wavelength range, leading to absolute MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? magnitude estimates consistent with a planetary-mass companion(2–5 𝑀 𝐽 considering the COND models and an age of 1–2 Myr;Baraffe et al. 2003).(vi) We fitted our spectrum with all available templates from theMontreal Spectral Library and found the best match with templatesof young mid-M dwarfs. This result is consistent with the measured 𝐽 -band spectral indices suggesting a spectral type of M5.5 ± 𝑅 𝐽 (considering all models within AIC-AIC min < 𝐾 𝛾 emission line, which would indicate on-going massaccretion.(viii) Our two leading hypotheses regarding the nature of thecompanion are: (1) an accreting protoplanet, from which we areprobing the fraction of the photosphere that is heated by accretionshocks; (2) an obscured stellar binary harbouring an edge-on discsuch that only a small fraction of its light, scattered from the discsurface, is reaching us.New observations at either shorter (including the H 𝛼 line)or longer wavelengths (mid-infrared and sub-mm) could confirmwhether the point source is an accreting substellar companion or anobscured stellar binary. ACKNOWLEDGMENTS
We thank Matthias Schreiber and Alejandro Melo for sharing thedata presented in Romero et al. (2012), Rebecca Jensen-Clem forsuggesting the use of a criterion based on the autocorrelation timefor MCMC convergence, and Julien Milli for useful discussions re-garding the degree of polarisation of companions. We acknowledgefunding from the Australian Research Council via DP180104235,FT130100034 and FT170100040. Part of this work has receivedfunding from the European Research Council (ERC) under the Eu-ropean Union’s Horizon 2020 research and innovation programme(grant agreement No 819155), and by the Wallonia-Brussels Fed-eration (grant for Concerted Research Actions. This project hasreceived funding from the European Union’s Horizon 2020 re-search and innovation programme under the Marie Skłodowska-Curie grant agreement No 823823 (DUSTBUSTERS). G-DM ac-knowledges the support of the German Science Foundation (DFG)priority program SPP 1992 “Exploring the Diversity of Extraso-lar Planets” (KU 2849/7-1) and from the Swiss National ScienceFoundation under grant BSSGI0_155816 “PlanetsInTime”. Parts ofthis work have been carried out within the framework of the NCCRPlanetS supported by the Swiss National Science Foundation. P.D.acknowledges the support of the French National Research Agencyin the framework of the Investissements d’Avenir program (ANR-15-IDEX-02), through the funding of the "Origin of Life" project ofthe Univ. Grenoble-Alpes. This work has made use of the SPHEREData Centre, jointly operated by OSUG/IPAG (Grenoble), PYTH-EAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observa-toire de Paris/LESIA (Paris), and Observatoire de Lyon/CRAL, andsupported by a grant from Labex OSUG@2020 (Investissementsd’avenir – ANR10 LABX56). This work has made use of datafrom the European Space Agency (ESA) mission
Gaia ( ), processed by the Gaia
Data Pro-cessing and Analysis Consortium (DPAC, ). Funding for the DPAChas been provided by national institutions, in particular the institu-tions participating in the
Gaia
Multilateral Agreement. This workhas made use of the Multi-modal Australian ScienceS Imaging andVisualisation Environment (MASSIVE) ( ).This research has benefitted from the Montreal Brown Dwarf andExoplanet Spectral Library, maintained by Jonathan Gagné.
DATA AVAILABILITY
The data underlying this article will be made available on CDS.
REFERENCES
Absil O., Mawet D., 2010, A&ARv, 18, 317Absil O., et al., 2013, A&A, 559, L12Akaike H., 1974, IEEE Transactions on Automatic Control, 19, 716Allard F., 2014, in Booth M., Matthews B. C., Graham J. R., eds, IAU Sym-posium Vol. 299, Exploring the Formation and Evolution of PlanetarySystems. pp 271–272, doi:10.1017/S1743921313008545Allard F., Homeier D., Freytag B., 2011, in Johns-Krull C., Browning M. K.,West A. A., eds, Astronomical Society of the Pacific Conference SeriesVol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems,and the Sun. p. 91 ( arXiv:1011.5405 )Allard F., Homeier D., Freytag B., 2012, Philosophical Transactions of theRoyal Society of London Series A, 370, 2765Allers K. N., et al., 2007, ApJ, 657, 511Amara A., Quanz S. P., 2012, MNRAS, 427, 948Aoyama Y., Marleau G.-D., Mordasini C., Ikoma M., 2020, arXiv e-prints,p. arXiv:2011.06608Auman Jason J., 1967, ApJS, 14, 171Baraffe I., Chabrier G., Barman T. S., Allard F., Hauschildt P. H., 2003,A&A, 402, 701Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577, A42Bate M. R., 2018, MNRAS, 475, 5618Beckwith S. V. W., Sargent A. I., Chini R. S., Guesten R., 1990, AJ, 99, 924Beuzit J.-L., et al., 2008, in Society of Photo-Optical Instrumentation Engi-neers (SPIE) Conference Series. , doi:10.1117/12.790120Bonnefoy M., et al., 2013, A&A, 555, A107Bonnefoy M., Chauvin G., Lagrange A.-M., Rojo P., Allard F., Pinte C.,Dumas C., Homeier D., 2014, A&A, 562, A127Boss A. P., 1998, ApJ, 503, 923Bouvier J., Dougados C., Alencar S. H. P., 2004, Ap&SS, 292, 659Bowler B. P., 2016, PASP, 128, 102001Brandt T. D., Dupuy T. J., Bowler B. P., 2019, AJ, 158, 140Brandt T. D., Dupuy T. J., Bowler B. P., Bardalez Gagliuffi D. C., FahertyJ., Brandt G. M., Michalik D., 2020, AJ, 160, 196Bresnahan D., et al., 2018, A&A, 615, A125Burnham K. P., Anderson D. R., eds, 2002, Information and LikelihoodTheory: A Basis for Model Selection and Inference. Springer New York,New York, NY, pp 49–97, doi:10.1007/978-0-387-22456-5_2, https://doi.org/10.1007/978-0-387-22456-5_2
Calcino J., Christiaens V., Price D. J., Pinte C., Davis T. M., van der MarelN., Cuello N., 2020, MNRAS, 498, 639Calvet N., Muzerolle J., Briceño C., Hernández J., Hartmann L., SaucedoJ. L., Gordon K. D., 2004, AJ, 128, 1294Canovas H., Rodenhuis M., Jeffers S. V., Min M., Keller C. U., 2011, A&A,531, A102Canovas H., et al., 2015, ApJ, 805, 21Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245Casassus S., 2016, PASA, 33, e013Cazzoletti P., et al., 2019, A&A, 626, A11Chauvin G., et al., 2012, A&A, 542, A41Chauvin G., et al., 2017, A&A, 605, L9Christiaens V., et al., 2018, A&A, 617, A37MNRAS , 1–22 (2021) V. Christiaens et al.
Christiaens V., et al., 2019a, MNRAS, 486, 5819Christiaens V., Cantalloube F., Casassus S., Price D. J., Absil O., Pinte C.,Girard J., Montesinos M., 2019b, ApJ, 877, L33Cieza L. A., et al., 2010, ApJ, 712, 925Claudi R. U., et al., 2008, in Ground-based and Airborne Instrumentationfor Astronomy II. p. 70143E, doi:10.1117/12.788366Cruz K. L., Reid I. N., 2002, AJ, 123, 2828Cuello N., et al., 2020, MNRAS, 491, 504Currie T., et al., 2019, ApJ, 877, L3Cutri R. M., et al., 2003, VizieR Online Data Catalog, p. II/246Delorme P., et al., 2017a, in Reylé C., Di Matteo P., Herpin F., Lagadec E.,Lançon A., Meliani Z., Royer F., eds, SF2A-2017: Proceedings of theAnnual meeting of the French Society of Astronomy and Astrophysics.p. Di ( arXiv:1712.06948 )Delorme P., et al., 2017b, A&A, 608, A79Dohlen K., et al., 2008, in McLean I. S., Casali M. M., eds, Society ofPhoto-Optical Instrumentation Engineers (SPIE) Conference Series Vol.7014, Ground-based and Airborne Instrumentation for Astronomy II. p.70143L, doi:10.1117/12.789786Dong R., et al., 2012, ApJ, 760, 111Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015, ApJ, 809, L5Dunham M. M., et al., 2015, ApJS, 220, 11Dzib S. A., Loinard L., Ortiz-León G. N., Rodríguez L. F., Galli P. A. B.,2018, ApJ, 867, 151Edwards S., Kwan J., Fischer W., Hillenbrand L., Finn K., Fedorenko K.,Feng W., 2013, ApJ, 778, 148Eriksson S. C., Asensio Torres R., Janson M., Aoyama Y., Marleau G.-D.,Bonnefoy M., Petrus S., 2020, A&A, 638, L6Espaillat C., et al., 2014, Protostars and Planets VI, pp 497–520Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125,306Gagné J., Lafrenière D., Doyon R., Malo L., Artigau É., 2014, ApJ, 783,121Gagné J., et al., 2015, ApJS, 219, 33Gaia Collaboration et al., 2016, A&A, 595, A1Gaia Collaboration et al., 2018, A&A, 616, A1Galicher R., et al., 2018, A&A, 615, A92Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457Ginski C., et al., 2018, A&A, 616, A79Girardi L., et al., 2012, Astrophysics and Space Science Proceedings, 26,165Gomez Gonzalez C. A., et al., 2017, AJ, 154, 7Goodman J., Weare J., 2010, Communications in Applied Mathematics andComputational Science, 5, 65Gravity Collaboration et al., 2019, A&A, 623, L11Greco J. P., Brandt T. D., 2016, ApJ, 833, 134Gressel O., Nelson R. P., Turner N. J., Ziegler U., 2013, ApJ, 779, 59Guidi G., et al., 2018, MNRAS, 479, 1505Gullbring E., Hartmann L., Briceño C., Calvet N., 1998, ApJ, 492, 323Gutermuth R. A., Megeath S. T., Myers P. C., Allen L. E., Pipher J. L., FazioG. G., 2009, ApJS, 184, 18Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., GirardJ. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749Haffert S. Y., et al., 2020, A&A, 640, L12Hartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135Hashimoto J., et al., 2011, ApJ, 729, L17Hashimoto J., Aoyama Y., Konishi M., Uyama T., Takasao S., Ikoma M.,Tanigawa T., 2020, arXiv e-prints, p. arXiv:2003.07922Hauschildt P. H., 1992, J. Quant. Spectrosc. Radiative Transfer, 47, 433Helling C., Woitke P., 2006, A&A, 455, 325Helling C., Dehn M., Woitke P., Hauschildt P. H., 2008, ApJ, 675, L105Herczeg G. J., Hillenbrand L. A., 2014, ApJ, 786, 97Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S., Pérez L.,2019, ApJ, 879, L25Johns-Krull C. M., 2007, ApJ, 664, 975Johns-Krull C. M., Valenti J. A., Koresko C., 1999, ApJ, 516, 900Keppler M., et al., 2018, A&A, 617, A44Keppler M., et al., 2020, A&A, 639, A62 Kratter K., Lodato G., 2016, ARA&A, 54, 271Kraus A. L., Ireland M. J., Cieza L. A., Hinkley S., Dupuy T. J., BowlerB. P., Liu M. C., 2014, ApJ, 781, 20Lagrange A.-M., et al., 2009, A&A, 493, L21Lagrange A.-M., et al., 2010, Science, 329, 57Launhardt R., et al., 2020, A&A, 635, A162Lenzen R., et al., 2003, in Iye M., Moorwood A. F. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841,Society of Photo-Optical Instrumentation Engineers (SPIE) ConferenceSeries. pp 944–952, doi:10.1117/12.460044Looper D. L., Bochanski J. J., Burgasser A. J., Mohanty S., Mamajek E. E.,Faherty J. K., West A. A., Pitts M. A., 2010, AJ, 140, 1486Lovelace R. V. E., Covey K. R., Lloyd J. P., 2011, AJ, 141, 51Luhman K. L., et al., 2008, ApJ, 675, 1375Maire A.-L., et al., 2016, in Ground-based and Airborne Instru-mentation for Astronomy VI. p. 990834 ( arXiv:1609.06681 ),doi:10.1117/12.2233013Marleau G. D., Cumming A., 2014, MNRAS, 437, 1378Marleau G.-D., Klahr H., Kuiper R., Mordasini C., 2017, ApJ, 836, 221Marleau G.-D., Coleman G. A. L., Leleu A., Mordasini C., 2019a, A&A,624, A20Marleau G.-D., Mordasini C., Kuiper R., 2019b, ApJ, 881, 144Marois C., Lafrenière D., Doyon R., Macintosh B., Nadeau D., 2006, ApJ,641, 556Marois C., Macintosh B., Barman T., Zuckerman B., Song I., Patience J.,Lafrenière D., Doyon R., 2008, Science, 322, 1348Marois C., Macintosh B., Véran J.-P., 2010, in Adaptive Optics Systems II.p. 77361J, doi:10.1117/12.857225Mawet D., et al., 2014, ApJ, 792, 97Mawet D., et al., 2016, in Adaptive Optics Systems V. p. 99090D,doi:10.1117/12.2233658McCabe C., Duchêne G., Pinte C., Stapelfeldt K. R., Ghez A. M., MénardF., 2011, ApJ, 727, 90Meyer M. R., Wilking B. A., 2009, PASP, 121, 350Milli J., Mouillet D., Lagrange A.-M., Boccaletti A., Mawet D., ChauvinG., Bonnefoy M., 2012, A&A, 545, A111Milli J., et al., 2017, A&A, 597, L2Mizuno H., 1980, Progress of Theoretical Physics, 64, 544Mordasini C., 2018, preprint, ( arXiv:1804.01532 )Mordasini C., Alibert Y., Klahr H., Henning T., 2012, A&A, 547, A111Mordasini C., Marleau G.-D., Mollière P., 2017, A&A, 608, A72Müller A., et al., 2018, A&A, 617, L2Muzerolle J., Hartmann L., Calvet N., 1998, AJ, 116, 2965Natta A., Testi L., Muzerolle J., Randich S., Comerón F., Persi P., 2004,A&A, 424, 603Naud M.-E., et al., 2014, ApJ, 787, 5Nelder & Mead 1965, Computer Journal, Vol 7, 308Neuhäuser R., et al., 2000, A&AS, 146, 323Nisini B., Antoniucci S., Giannini T., Lorenzetti D., 2005, A&A, 429, 543Nowak M., et al., 2020, A&A, 642, L2Nutter D. J., Ward-Thompson D., André P., 2005, MNRAS, 357, 975Olofsson J., et al., 2013, A&A, 552, A4Owen J. E., 2016, Publ. Astron. Soc. Australia, 33, e005Paardekooper S.-J., Johansen A., 2018, Space Sci. Rev., 214, 38Pairet B., Cantalloube F., Gomez Gonzalez C. A., Absil O., Jacques L.,2019, MNRAS, 487, 2262Pairet B., Cantalloube F., Jacques L., 2020, arXiv e-prints, p.arXiv:2008.05170Peterson D. E., et al., 2011, ApJS, 194, 43Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M.,Greenzweig Y., 1996, Icarus, 124, 62Price D. J., et al., 2018, MNRAS, 477, 1270Quanz S. P., Amara A., Meyer M. R., Kenworthy M. A., Kasper M., GirardJ. H., 2013, ApJ, 766, L1Rich E. A., et al., 2019, ApJ, 875, 38Romero G. A., Schreiber M. R., Cieza L. A., Rebassa-Mansergas A., MerínB., Smith Castelli A. V., Allen L. E., Morrell N., 2012, ApJ, 749, 79Rousset G., et al., 2003, in Wizinowich P. L., Bonaccini D., eds, SocietyMNRAS000
Christiaens V., et al., 2019a, MNRAS, 486, 5819Christiaens V., Cantalloube F., Casassus S., Price D. J., Absil O., Pinte C.,Girard J., Montesinos M., 2019b, ApJ, 877, L33Cieza L. A., et al., 2010, ApJ, 712, 925Claudi R. U., et al., 2008, in Ground-based and Airborne Instrumentationfor Astronomy II. p. 70143E, doi:10.1117/12.788366Cruz K. L., Reid I. N., 2002, AJ, 123, 2828Cuello N., et al., 2020, MNRAS, 491, 504Currie T., et al., 2019, ApJ, 877, L3Cutri R. M., et al., 2003, VizieR Online Data Catalog, p. II/246Delorme P., et al., 2017a, in Reylé C., Di Matteo P., Herpin F., Lagadec E.,Lançon A., Meliani Z., Royer F., eds, SF2A-2017: Proceedings of theAnnual meeting of the French Society of Astronomy and Astrophysics.p. Di ( arXiv:1712.06948 )Delorme P., et al., 2017b, A&A, 608, A79Dohlen K., et al., 2008, in McLean I. S., Casali M. M., eds, Society ofPhoto-Optical Instrumentation Engineers (SPIE) Conference Series Vol.7014, Ground-based and Airborne Instrumentation for Astronomy II. p.70143L, doi:10.1117/12.789786Dong R., et al., 2012, ApJ, 760, 111Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015, ApJ, 809, L5Dunham M. M., et al., 2015, ApJS, 220, 11Dzib S. A., Loinard L., Ortiz-León G. N., Rodríguez L. F., Galli P. A. B.,2018, ApJ, 867, 151Edwards S., Kwan J., Fischer W., Hillenbrand L., Finn K., Fedorenko K.,Feng W., 2013, ApJ, 778, 148Eriksson S. C., Asensio Torres R., Janson M., Aoyama Y., Marleau G.-D.,Bonnefoy M., Petrus S., 2020, A&A, 638, L6Espaillat C., et al., 2014, Protostars and Planets VI, pp 497–520Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125,306Gagné J., Lafrenière D., Doyon R., Malo L., Artigau É., 2014, ApJ, 783,121Gagné J., et al., 2015, ApJS, 219, 33Gaia Collaboration et al., 2016, A&A, 595, A1Gaia Collaboration et al., 2018, A&A, 616, A1Galicher R., et al., 2018, A&A, 615, A92Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457Ginski C., et al., 2018, A&A, 616, A79Girardi L., et al., 2012, Astrophysics and Space Science Proceedings, 26,165Gomez Gonzalez C. A., et al., 2017, AJ, 154, 7Goodman J., Weare J., 2010, Communications in Applied Mathematics andComputational Science, 5, 65Gravity Collaboration et al., 2019, A&A, 623, L11Greco J. P., Brandt T. D., 2016, ApJ, 833, 134Gressel O., Nelson R. P., Turner N. J., Ziegler U., 2013, ApJ, 779, 59Guidi G., et al., 2018, MNRAS, 479, 1505Gullbring E., Hartmann L., Briceño C., Calvet N., 1998, ApJ, 492, 323Gutermuth R. A., Megeath S. T., Myers P. C., Allen L. E., Pipher J. L., FazioG. G., 2009, ApJS, 184, 18Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., GirardJ. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749Haffert S. Y., et al., 2020, A&A, 640, L12Hartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135Hashimoto J., et al., 2011, ApJ, 729, L17Hashimoto J., Aoyama Y., Konishi M., Uyama T., Takasao S., Ikoma M.,Tanigawa T., 2020, arXiv e-prints, p. arXiv:2003.07922Hauschildt P. H., 1992, J. Quant. Spectrosc. Radiative Transfer, 47, 433Helling C., Woitke P., 2006, A&A, 455, 325Helling C., Dehn M., Woitke P., Hauschildt P. H., 2008, ApJ, 675, L105Herczeg G. J., Hillenbrand L. A., 2014, ApJ, 786, 97Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S., Pérez L.,2019, ApJ, 879, L25Johns-Krull C. M., 2007, ApJ, 664, 975Johns-Krull C. M., Valenti J. A., Koresko C., 1999, ApJ, 516, 900Keppler M., et al., 2018, A&A, 617, A44Keppler M., et al., 2020, A&A, 639, A62 Kratter K., Lodato G., 2016, ARA&A, 54, 271Kraus A. L., Ireland M. J., Cieza L. A., Hinkley S., Dupuy T. J., BowlerB. P., Liu M. C., 2014, ApJ, 781, 20Lagrange A.-M., et al., 2009, A&A, 493, L21Lagrange A.-M., et al., 2010, Science, 329, 57Launhardt R., et al., 2020, A&A, 635, A162Lenzen R., et al., 2003, in Iye M., Moorwood A. F. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841,Society of Photo-Optical Instrumentation Engineers (SPIE) ConferenceSeries. pp 944–952, doi:10.1117/12.460044Looper D. L., Bochanski J. J., Burgasser A. J., Mohanty S., Mamajek E. E.,Faherty J. K., West A. A., Pitts M. A., 2010, AJ, 140, 1486Lovelace R. V. E., Covey K. R., Lloyd J. P., 2011, AJ, 141, 51Luhman K. L., et al., 2008, ApJ, 675, 1375Maire A.-L., et al., 2016, in Ground-based and Airborne Instru-mentation for Astronomy VI. p. 990834 ( arXiv:1609.06681 ),doi:10.1117/12.2233013Marleau G. D., Cumming A., 2014, MNRAS, 437, 1378Marleau G.-D., Klahr H., Kuiper R., Mordasini C., 2017, ApJ, 836, 221Marleau G.-D., Coleman G. A. L., Leleu A., Mordasini C., 2019a, A&A,624, A20Marleau G.-D., Mordasini C., Kuiper R., 2019b, ApJ, 881, 144Marois C., Lafrenière D., Doyon R., Macintosh B., Nadeau D., 2006, ApJ,641, 556Marois C., Macintosh B., Barman T., Zuckerman B., Song I., Patience J.,Lafrenière D., Doyon R., 2008, Science, 322, 1348Marois C., Macintosh B., Véran J.-P., 2010, in Adaptive Optics Systems II.p. 77361J, doi:10.1117/12.857225Mawet D., et al., 2014, ApJ, 792, 97Mawet D., et al., 2016, in Adaptive Optics Systems V. p. 99090D,doi:10.1117/12.2233658McCabe C., Duchêne G., Pinte C., Stapelfeldt K. R., Ghez A. M., MénardF., 2011, ApJ, 727, 90Meyer M. R., Wilking B. A., 2009, PASP, 121, 350Milli J., Mouillet D., Lagrange A.-M., Boccaletti A., Mawet D., ChauvinG., Bonnefoy M., 2012, A&A, 545, A111Milli J., et al., 2017, A&A, 597, L2Mizuno H., 1980, Progress of Theoretical Physics, 64, 544Mordasini C., 2018, preprint, ( arXiv:1804.01532 )Mordasini C., Alibert Y., Klahr H., Henning T., 2012, A&A, 547, A111Mordasini C., Marleau G.-D., Mollière P., 2017, A&A, 608, A72Müller A., et al., 2018, A&A, 617, L2Muzerolle J., Hartmann L., Calvet N., 1998, AJ, 116, 2965Natta A., Testi L., Muzerolle J., Randich S., Comerón F., Persi P., 2004,A&A, 424, 603Naud M.-E., et al., 2014, ApJ, 787, 5Nelder & Mead 1965, Computer Journal, Vol 7, 308Neuhäuser R., et al., 2000, A&AS, 146, 323Nisini B., Antoniucci S., Giannini T., Lorenzetti D., 2005, A&A, 429, 543Nowak M., et al., 2020, A&A, 642, L2Nutter D. J., Ward-Thompson D., André P., 2005, MNRAS, 357, 975Olofsson J., et al., 2013, A&A, 552, A4Owen J. E., 2016, Publ. Astron. Soc. Australia, 33, e005Paardekooper S.-J., Johansen A., 2018, Space Sci. Rev., 214, 38Pairet B., Cantalloube F., Gomez Gonzalez C. A., Absil O., Jacques L.,2019, MNRAS, 487, 2262Pairet B., Cantalloube F., Jacques L., 2020, arXiv e-prints, p.arXiv:2008.05170Peterson D. E., et al., 2011, ApJS, 194, 43Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M.,Greenzweig Y., 1996, Icarus, 124, 62Price D. J., et al., 2018, MNRAS, 477, 1270Quanz S. P., Amara A., Meyer M. R., Kenworthy M. A., Kasper M., GirardJ. H., 2013, ApJ, 766, L1Rich E. A., et al., 2019, ApJ, 875, 38Romero G. A., Schreiber M. R., Cieza L. A., Rebassa-Mansergas A., MerínB., Smith Castelli A. V., Allen L. E., Morrell N., 2012, ApJ, 749, 79Rousset G., et al., 2003, in Wizinowich P. L., Bonaccini D., eds, SocietyMNRAS000 , 1–22 (2021) rA-9 b, B or BKG? of Photo-Optical Instrumentation Engineers (SPIE) Conference SeriesVol. 4839, Society of Photo-Optical Instrumentation Engineers (SPIE)Conference Series. pp 140–149, doi:10.1117/12.459332Sallum S., et al., 2015, Nature, 527, 342Schmid H. M., Joos F., Tschan D., 2006, A&A, 452, 657Sicilia-Aguilar A., Henning T., Juhász A., Bouwman J., Garmire G., GarmireA., 2008, ApJ, 687, 1145Sicilia-Aguilar A., Henning T., Kainulainen J., Roccatagliata V., 2011, ApJ,736, 137Sicilia-Aguilar A., Henning T., Linz H., André P., Stutz A., Eiroa C., WhiteG. J., 2013, A&A, 551, A34Sokal K. R., Deen C. P., Mace G. N., Lee J.-J., Oh H., Kim H., Kidder B. T.,Jaffe D. T., 2018, ApJ, 853, 120Soummer R., Pueyo L., Larkin J., 2012, ApJ, 755, L28Sparks W., Ford H., 2002, ApJ, 578, 543Spiegel D. S., Burrows A., 2012, ApJ, 745, 174Spina L., et al., 2017, A&A, 601, A70Stapelfeldt K. R., Krist J. E., Ménard F., Bouvier J., Padgett D. L., BurrowsC. J., 1998, ApJ, 502, L65Stolker T., et al., 2020a, A&A, 635, A182Stolker T., Marleau G. D., Cugno G., Mollière P., Quanz S. P., TodorovK. O., Kühn J., 2020b, A&A, 644, A13Szulágyi J., 2017, ApJ, 842, 103Taylor K. N. R., Storey J. W. V., 1984, MNRAS, 209, 5PTognelli E., Prada Moroni P. G., Degl’Innocenti S., 2011, A&A, 533, A109Ubeira-Gabellini M. G., Christiaens V., Lodato G., Ancker M. v. d., FedeleD., Manara C. F., Price D. J., 2020, ApJ, 890, L8Wang H., Mundt R., Henning T., Apai D., 2004, ApJ, 617, 1191Wang J. J., et al., 2020, AJ, 159, 263Wang J. J., et al., 2021, arXiv e-prints, p. arXiv:2101.04187Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296Wertz O., Absil O., Gómez González C. A., Milli J., Girard J. H., MawetD., Pueyo L., 2017, A&A, 598, A83Williams J. P., Cieza L. A., 2011, ARA&A, 49, 67Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409Woitke P., Helling C., 2003, A&A, 399, 297Woitke P., Helling C., 2004, A&A, 414, 335Wright E. L., et al., 2010, AJ, 140, 1868Wu Y.-L., Sheehan P. D., 2017, ApJ, 846, L26Zhu Z., 2015, ApJ, 799, 16Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88van der Bliek N. S., Manfroid J., Bouchet P., 1996, A&AS, 119, 547van der Marel N., Cazzoletti P., Pinilla P., Garufi A., 2016, ApJ, 832, 178van der Marel N., et al., 2021, AJ, 161, 33 APPENDIX A: VLT/NACO REDUCTION PIPELINE FORNON-CORONAGRAPHIC DATA (i) First the pipeline automatically finds the stellar position ineach science cube. This is done by a) subtracting the median of theclosest three cubes in time (excluding the cube in question), whichsubtracts an estimate of both DARK current and sky background; b)removing spatial frequencies corresponding to either bad pixels orlarge scale detector variation (i.e. using a low-pass and a high-passfilter); c) looking for maxima in the median image of each cube. Werecord the approximate ( 𝑥 , 𝑦 ) position of the star on the detector andthe quadrant in which it is located.(ii) We then subtract the sky from the original frames of eachcube, using the principal component analysis (PCA) based sky sub-traction algorithm implemented in vip. For each science cube, thePCA library was set to the medians of each individual cube wherethe star is located in a different quadrant from the one considered.Both science and sky cubes are cropped to the quadrant of the de-tector in which the star is located. We tested 1 to 30 subtractedprincipal components ( 𝑛 pc ), and noted no further improvement in dark current and sky background residuals for 𝑛 pc larger than 5. Wehence adopted the latter value.(iii) The pipeline then calculates a master flat field from raw flatsacquired during twilight at three different telescope altitudes, andapply it to all science cubes. It also computes a static bad pixel mapbased on flat values smaller than 0.67 and larger than 1.5.(iv) We correct for NaN values, and subsequently bad pixels,using an iterative sigma filter algorithm implemented in vip anddesigned to correct for clusters of bad pixels. The first pass of badpixel correction uses the static bad pixel map inferred from the flatfield, while the second pass looks iteratively for 8- 𝜎 outliers tocorrect for cosmic rays. Bad pixels are replaced by the median ofgood neighbouring pixels.(v) We subsequently find the centroid of the star by fitting a 2DMoffat function and shift each frame so that the star would fallexactly on the central pixel of each image (images are cropped toodd dimensions).(vi) All recentered images are then gathered into a single mastercube. For each frame of the master cube, the associated derotationangle required to align North up and East left is found by interpo-lation of the parallactic angle at start and end of the original cube itpertains to.(vii) From which we identify and remove bad frames, either cor-responding to the opening/decrease in performance of the adaptiveoptics loop, or jittering during the integration which elongated thePSF. Bad frame trimming is based on the cross-correlation of eachframe to the median of all frames in the cube. For each frame, thethreshold used to consider a frame bad is a computed Pearson cor-relation factor (with the median frame) lower than 0.8, as calculatedin a 7 FWHM x 7 FWHM cropped window centered on the star. Forboth NACO datasets, this removed ∼
5% of all images.(viii) Given the large number of images (over 26,000 and 55,000for the 2017 and 2018 datasets, respectively), we then median-combine consecutive frames together 10 by 10 and 16 by 16 for the2017 and 2018 datasets, respectively.(ix) The FWHM of the stellar point-spread function is then es-timated by fitting a 2D gaussian profile to the median image of thewhole cube. The flux of the star is then calculated by integrationover a 1-FWHM aperture.(x) Finally, we post-process the calibrated cubes using bothmedian-ADI (Marois et al. 2006) in 5 (cid:48)(cid:48) × (cid:48)(cid:48) frames and PCA-ADI(e.g. Amara & Quanz 2012) in 2 (cid:48)(cid:48) × (cid:48)(cid:48) frames, as implemented invip. For PCA-ADI, a range of 1 to 100 𝑛 pc was explored. APPENDIX B: IRDIS REDUCTION PIPELINE FORNON-CORONAGRAPHIC DATA (i) The pipeline first calculates master sky background imagesusing esorex’s sph_ird_sky_bg recipe, upon provision of the rawsky images obtained during the sequence.(ii) A master flat field and static bad pixel map are then calculatedwith the sph_ird_instrument_flat recipe, using raw flats andcorresponding raw darks.(iii) The master sky is subtracted and the flat field is divided fromall good pixels of science cubes using the sph_ird_science_dbi recipe. Since this does not systematically yield an average back-ground level of zero, we complemented by an additional manualsubtraction of the residual sky so that the median pixel intensitiesat > (cid:48)(cid:48) from the star is null.(iv) Next, the pipeline uses vip routines to correct for bad pixelsusing an iterative sigma filter algorithm. A first pass corrects for the MNRAS , 1–22 (2021) V. Christiaens et al. static bad pixels identified in the master flat field, and a second passcorrects for all residual 8 𝜎 outliers.(v) A 2D Moffat profile is subsequently fit to the stellar PSF ineach frame in order to find the centroid, and all images are shiftedfor the star to fall exactly on the central pixel of odd-size frames.(vi) All centered cubes are then collated in a single master cube,and corresponding derotation angles to align north up and east leftare calculated. The derotation angles are interpolated for each frameof each cube, based on the parallactic angles at start and end of eachcube, and consider the true north value of -1.75 ± ◦ measured inMaire et al. (2016).(vii) Bad frames are then identified and removed according toboth pixel intensities (rejecting stellar fluxes 1 𝜎 below the medianflux) and Pearson correlation coefficient calculated with respect tothe median frame.(viii) The anamorphism measured in Maire et al. (2016) is thencorrected by rescaling the image along the 𝑦 dimension using afourth-order Lanczos interpolation.(ix) The pipeline finally post-processes the calibrated cube usingmedian-ADI on full frames, and both PCA-ADI and sPCA (Absilet al. 2013) on 2 (cid:48)(cid:48) × (cid:48)(cid:48) cropped frames. sPCA performs PCA-ADIin concentric annuli, with the PCA library for the annulus of eachimage built from the same annulus in other images of the cube wherea putative planet would have rotated by more than a given threshold.The PCA-ADI algorithms are run with 1 to 50 𝑛 pc subtracted, andan angular threshold ranging from 0.5 × FWHM (for the innermostannulus) to 1 × FWHM (for the outermost annulus) linear motionwas chosen for sPCA.
APPENDIX C: IFS REDUCTION PIPELINE FORNON-CORONAGRAPHIC DATA (i) The pipeline first calculates master darks for all images (sci-ence or other calibrations) with different integration times usingthe sph_ifs_master_dark recipe. This is necessary because theesorex IFS recipes (version 3.13.2) do not appear to scale masterdarks when applied to images obtained with different integrationtimes. Master darks are subtracted manually to the raw flats, pairingthem based on integration time.(ii) Master detector flats are then calculated using the sph_ifs_master_detector_flat recipe. This is done in foursteps, where the output of each step is used as input in the followingstep: first a preamplifier flat is calculated using broad-band lampraw flats; second large-scale coloured flats are calculated for eachof the 4 narrow-band flat lamps; third a large-scale white lamp flatis calculated; and fourth a small-scale coloured flat is calculated af-ter removing the large scale structures (captured with a smoothinglength of 5 pixels). In all subsequent recipes, we provided eitherlarge-scale coloured flats and/or the small-scale white flat as detec-tor flat inputs for optimal performance.(iii) Next, master sky background images are calculated from thesky cubes acquired during the observation, using the recipe.(iv) The spectra positions on the IFS detector are then determinedusing the sph_ifs_spectra_positions recipe. We changed thedefault value for the distortion option and set it to False, as wenoticed letting it to True led to significant negative/positive parallelstripes in some of the spectral channels of the final cubes.(v) Both the master detector flats and spectra positions are sub-sequently provided as input to the sph_ifs_instrument_flat recipe in order to compute a total instrument flat.(vi) The exact wavelength solution is then found by providing both the outputs of the above steps and raw wave calibration filesusing the sph_ifs_wave_calib recipe.(vii) A master IFU flat is calculated by passing the outputs ofprevious steps to the sph_ifs_instrument_flat recipe.(viii) All science cubes are then calibrated by providing the out-puts of steps ii to vii to the sph_ifs_science_dr recipe. Sincesky subtraction is sub-optimal, we manually correct for the residualsky level in order to have a median level of zero in 0 . (cid:48)(cid:48) 𝑦 and 𝑥 owing to the rotation of the field-of-view.(x) Finally, the pipeline post-processes the calibrated 4D IFScube leveraging on spectral differential imaging (SDI; Sparks &Ford 2002) and angular differential imaging. More specifically: (a)PCA-SDI was applied on each individual spectral cube, and the SDIimages are then derotated and median-combined (simply referredto as PCA-SDI throughout this paper); (b) PCA-ADI was appliedat each wavelength on the 3D cubes composed of each individualspectral channel sampled in the temporal direction; (c) and PCA-ASDI, combining both the radial and azimuthal diversity of SDIand ADI together, was applied to reach the highest contrast at theexpense of possible self-subtraction of extended disc features. Asopposed to the version of PCA-ASDI used in Christiaens et al.(2019a) on SINFONI data the algorithm only builds a single PCAlibrary containing both spectral and angular diversity. Our tests onthese SPHERE/IFS data suggest that PCA-ASDI in two steps reducethe SNR of the companion and appear to significantly lower thealgorithmic throughput without gain in achieved contrast. A rangeof 1 to 10 𝑛 pc was explored for both PCA-SDI and PCA-ASDI. APPENDIX D: IMPROVEMENTS TO MCMC-NEGFC
In order for the MCMC algorithm to work (and converge) properly,we now use the following log-probability:log L( 𝐷 | 𝑀 ) = − ∑︁ 𝑖 ( 𝐼 𝑖 − 𝜇 ) 𝜎 (D1)where 𝑖 are the indices of all the pixels in a 1.5-FWHM radiusaperture centered on the location of our initial guess on the po-sition of the companion; 𝐼 𝑖 are the pixel intensities measured inthat aperture in the post-processed PCA-ADI image; 𝜇 and 𝜎 arethe mean intensity and standard deviation measured in a 3-FWHMwide annulus at the same radial separation as the companion butexcluding a region in azimuth near the location of the companion.The exclusion region is set to [PA 𝑏 - Δ PA, PA 𝑏 + Δ PA], where PA 𝑏 isour initial estimate on the position angle of the companion, and Δ PAis the range of parallactic angles covered by the ADI sequence (seeTable 2). We argue that Equation D1 leads to a more robust estimateof the companion flux than the original expression (missing 𝜇 and 𝜎 ). If the speckle residuals are not perfectly subtracted, as can bethe case for a small 𝑛 pc subtracted, variable adaptive optics perfor-mance throughout the sequence, and/or companions very close tothe central star, the mean residual level 𝜇 at the separation of thecompanion can be non-null (e.g. Christiaens et al. 2018). It is thusnecessary to subtract that component. Scaling by 𝜎 measured in anannulus that avoids the area around the point source also allows toinclude the speckle noise as an additional source of variance in theposterior distributions. Furthermore, for very faint companions in MNRAS000
In order for the MCMC algorithm to work (and converge) properly,we now use the following log-probability:log L( 𝐷 | 𝑀 ) = − ∑︁ 𝑖 ( 𝐼 𝑖 − 𝜇 ) 𝜎 (D1)where 𝑖 are the indices of all the pixels in a 1.5-FWHM radiusaperture centered on the location of our initial guess on the po-sition of the companion; 𝐼 𝑖 are the pixel intensities measured inthat aperture in the post-processed PCA-ADI image; 𝜇 and 𝜎 arethe mean intensity and standard deviation measured in a 3-FWHMwide annulus at the same radial separation as the companion butexcluding a region in azimuth near the location of the companion.The exclusion region is set to [PA 𝑏 - Δ PA, PA 𝑏 + Δ PA], where PA 𝑏 isour initial estimate on the position angle of the companion, and Δ PAis the range of parallactic angles covered by the ADI sequence (seeTable 2). We argue that Equation D1 leads to a more robust estimateof the companion flux than the original expression (missing 𝜇 and 𝜎 ). If the speckle residuals are not perfectly subtracted, as can bethe case for a small 𝑛 pc subtracted, variable adaptive optics perfor-mance throughout the sequence, and/or companions very close tothe central star, the mean residual level 𝜇 at the separation of thecompanion can be non-null (e.g. Christiaens et al. 2018). It is thusnecessary to subtract that component. Scaling by 𝜎 measured in anannulus that avoids the area around the point source also allows toinclude the speckle noise as an additional source of variance in theposterior distributions. Furthermore, for very faint companions in MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? NACO 2017 IFS 2019 (channel 1) IRDIS 2019 (K1) plate scale: 27.1 1.1 mas/pxtrue north: -0.58 0.30º (2017)true north: -0.58 1.00º (2018) ±±± plate scale: 7.46 0.02 mas/pxtrue north: -102.18 0.13ºpupil offset: -135.99 0.11º ± ±± plate scale: 12.267 0.010 mas/pxtrue north: -1.75 0.08º ±± Figure D1.
Corner plots for the radial separation, azimuth angle and contrast of the companion candidate retrieved by MCMC-NEGFC on the NACO 2017,IFS 2019 (first channel; lowest SNR) and IRDIS 2019 ( 𝐾 𝑟 are in pixels and the azimuthal angles 𝜃 aremeasured counter-clockwise from the positive 𝑥 axis. We also show the systematic uncertainties considered for the calculation of final astrometric uncertaintiesfor each instrument (NACO: Milli et al. 2017; SPHERE: Maire et al. 2016). raw detector units (e.g. in ADI datacubes that have been normalisedby the integration time), the absence of scaling by 𝜎 would leadto large log L( 𝐷 | 𝑀 ) likelihoods whether the companion is wellsubtracted or not subtracted at all. In turn this leads to an unreliableposterior distribution for the flux, hence associated best estimateand uncertainty.The MCMC-NEGFC algorithm stops once a given conver-gence criterion is met. In the original version of the algorithm, theconvergence criterion relied on a Gelman-Rubin test; i.e. based ona comparison of the variance calculated from two sections of thechains (Gelman & Rubin 1992). However, since the samples in aMarkov chain are not independent, the Gelman-Rubin test is inade-quate and may break the chain too early, resulting in underestimatedvariances of the posterior distributions, hence uncertainties. There-fore, we now use a test based on the integrated auto-correlationtime 𝜏 (e.g. Goodman & Weare 2010; Foreman-Mackey et al.2013). The estimate of 𝜏 is performed at regular intervals as thechain progresses, and becomes reliable if the number of samples 𝑁 on which it is computed is sufficiently large compared to theestimated 𝜏 ( 𝑁 >> 𝜏 ). In practice, we considered a threshold of 𝑁 / max ( 𝜏 𝑓 ) >
50 where 𝜏 𝑓 is the autocorrelation time for parame-ter 𝑓 (i.e. 𝑟 , PA or contrast), as recommended in the documentationof emcee (Foreman-Mackey et al. 2013). Since the sampling uncer-tainty on the true variance of a parameter is ∝ √︁ 𝜏 𝑓 / 𝑁 , the criterionwe chose also implies a relative accuracy of (cid:46)
14% on the variances(and thereby uncertainties) inferred for each parameter. Using thisconvergence criterion, we noticed that a significantly larger num-ber of steps was required before the convergence criterion was metcompared to the Gelman-Rubin-based criterion: between 1500 and2500 iterations were required for all our datasets, while the Gelman-Rubin-based criterion typically suggested convergence within 200steps. We used 128 walkers for all datasets and a burn-in factor of0.5, leading to a total of > 𝑟 , PA and contrast of the point source in each dataset.Our third modification is motivated by the varying quality ofthe adaptive optics correction throughout all our ADI sequences,which resulted in significant scatter for the measured stellar fluxes(a proxy for the Strehl ratio). Only the NACO 2017 dataset showed a low relative standard deviation of 3.2% for stellar fluxes measuredin all frames (with respect to the median flux). The NACO 2018,IFS 2019 and IRDIS 2019 (unsaturated) datasets showed relativestandard deviations of 24.7%, 13.9% and 8.7% (after bad frameremoval), with individual flux measurements varying by up to afactor 2. Fortunately, since the point source was recovered in allour unsaturated datasets, we could measure the stellar flux in eachindividual image, and use that information to inject the flux of thenegative companion in each image proportionally. The injected fluxin frame 𝑖 is: 𝐹 𝑖 = 𝐹 × 𝐹 ∗ ,𝑖 𝐹 ∗ , med (D2)where 𝐹 is the companion candidate flux sampled by the MCMC-NEGFC algorithm, 𝐹 ∗ ,𝑖 is the stellar flux measured in frame 𝑖 and 𝐹 ∗ , med is the median stellar flux. This modification allowed us toget up to an order of magnitude improvement in accuracy for theestimated contrast of the point source with respect to the star inthe different datasets. Since for coronagraphic datasets one doesnot have the luxury of knowing how the stellar flux (hence thecompanion flux) varies throughout the observing sequence, thismodification was only implemented as an additional option in vip’smcmc_negfc_sampling function.Figure D1 shows three examples of the results obtained byMCMC-NEGFC among the 43 ADI sequences considered – onefor each filter or spectral channel of the NACO, IFS and IRDISdatasets. We selected the best NACO 𝐿 (cid:48) dataset (2017), the firstchannel of the IFS dataset (lowest SNR for the companion), and the 𝐾 APPENDIX E: IMAGES OBTAINED WITH PCA-SDI,PCA-ASDI AND MAYO
Figure E1 shows the post-processed IFS images obtained with PCA-SDI and PCA-ASDI for increasing number of principal componentssubtracted. The purpose of the figure is to evaluate the reliabilityof the spiral pattern. As expected for an authentic disc feature, the
MNRAS , 1–22 (2021) V. Christiaens et al. a r c s ec
50 auPCA-SDI (npc = 1) 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8arcsec-0.8-0.6-0.4-0.20.00.20.40.60.8 a r c s ec
50 auPCA-SDI (npc = 2) 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8arcsec-0.8-0.6-0.4-0.20.00.20.40.60.8 a r c s ec
50 auPCA-SDI (npc = 3) 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8arcsec-0.8-0.6-0.4-0.20.00.20.40.60.8 a r c s ec
50 auPCA-SDI (npc = 4) a r c s ec
50 auPCA-ASDI (npc = 2) 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8arcsec-0.8-0.6-0.4-0.20.00.20.40.60.8 a r c s ec
50 auPCA-ASDI (npc = 3) 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8arcsec-0.8-0.6-0.4-0.20.00.20.40.60.8 a r c s ec
50 auPCA-ASDI (npc = 4) 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8arcsec-0.8-0.6-0.4-0.20.00.20.40.60.8 a r c s ec
50 auPCA-ASDI (npc = 5)
Figure E1.
PCA-SDI (top row) and PCA-ASDI (bottom row) images obtained on the SPHERE/IFS dataset. A tentative spiral pattern is detected, with a possibleprimary arm pointing towards the point source. a r c s ec
50 auMAYO (K1+K2)
SPHERE/IRDIS K1K2 (2019-09-30)
Figure E2.
IRDIS 𝐾 𝐾 morphology of the signal is preserved but gets progressively self-subtracted for increasing 𝑛 pc . Compared to the PCA-SDI images,the PCA-ASDI images removed all the azimuthally extended signalsdue to angular differential imaging. Since we were only interestedin testing the reliability of the spiral pattern, we ran the PCA-ASDIalgorithm using the crop_ifs option of vip’s pca routine, whichcrops frames after rescaling in order to both save memory anddecrease computation time. Cropping removes different amount ofinformation radially in different spectral channels and accounts forthe abrupt edges of the cropped field in the final median-combinedframe, however it does not affect the signal at shorter separation. We also applied mayonnaise on the IRDIS dataset (Pairet et al.2020). Contrary to PCA, mayonnaise projects sparse (planet-like)and extended (disc-like) signals on different bases, which allowsboth components to be restored and disentangled. mayonnaise re-covered sparse emission from the companion, and possibly extendedemission from the protoplanetary disc (Figure E2). However, giventhe absence of spirals in the IRDIS image and their low SNR in theIFS image ( ∼ 𝑛 pc = APPENDIX F: CORNER PLOT OF CRA-9
The corner plot showing the most likely physical parameters for thestar inferred by specfit is presented in Figure F1.
APPENDIX G:
SPECFIT
TESTS ON THE COMPANIONSPECTRUM
Figure G1 shows the best-fit model obtained after running specfitusing a gaussian prior of 𝜇 = . 𝜎 = . ( 𝑔 ) withthe BT-SETTL grid of models. The favoured parameters are 𝑇 𝑒 = ( 𝑔 ) = . 𝑅 𝑏 = . 𝑅 𝐽 and 𝐴 𝑉 = . Δ AIC ∼ 𝐻 band measurements being slightlyoverpredicted.The BT-SETTL and DRIFT-PHOENIX models favoured byspecfit when the photometric radius is fixed to 1.8 𝑅 𝐽 are shownin Figure G2. We did not test BT-DUSTY models with such con-straint due to the incompleteness of this grid at lower temperaturesthan 3000 K. The plots show that this constraint leads to visuallypoor fits. This is also conveyed by the Δ AIC values of ∼ ∼ MNRAS000
Figure G1 shows the best-fit model obtained after running specfitusing a gaussian prior of 𝜇 = . 𝜎 = . ( 𝑔 ) withthe BT-SETTL grid of models. The favoured parameters are 𝑇 𝑒 = ( 𝑔 ) = . 𝑅 𝑏 = . 𝑅 𝐽 and 𝐴 𝑉 = . Δ AIC ∼ 𝐻 band measurements being slightlyoverpredicted.The BT-SETTL and DRIFT-PHOENIX models favoured byspecfit when the photometric radius is fixed to 1.8 𝑅 𝐽 are shownin Figure G2. We did not test BT-DUSTY models with such con-straint due to the incompleteness of this grid at lower temperaturesthan 3000 K. The plots show that this constraint leads to visuallypoor fits. This is also conveyed by the Δ AIC values of ∼ ∼ MNRAS000 , 1–22 (2021) rA-9 b, B or BKG? T e = +189−137 K . . . . l og ( g ) log( g ) = +0.2−0.2 . . . . R * R * = +0.09−0.14 R ⊙ . . . . A V A V = +0.3−0.2 mag T e . . . . R V . . . . log( g ) . . . . R * . . . . A V . . . . R V R V = +1.7−0.5 Figure F1.
Corner plot showing the posterior distribution for the parametersderived by specfit for CrA-9. We note minor degeneracies (i) between theeffective temperature 𝑇 𝑒 , the radius 𝑅 ∗ and the optical extinction 𝐴 𝑉 , and(ii) between 𝐴 𝑉 and the total-to-selective optical extinction ratio 𝑅 𝑉 . respectively. For the BT-SETTL fit, we notice that specfit con-verged on three (equally bad) clusters of parametric solutions witheffective temperature lower than 2000 K, while a single kind ofsolution at an effective temperature of ∼ 𝑇 𝑒 = ( 𝑔 ) = . 𝑅 𝑏 = . 𝑅 𝐽 and 𝐴 𝑉 = . 𝑇 𝑒 = ( 𝑔 ) = . 𝑅 𝑏 = . 𝑅 𝐽 and 𝐴 𝑉 = . This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS , 1–22 (2021) V. Christiaens et al. μ m)0.00.51.01.52.02.5 λ F λ ( W m − μ m − ) T e = 3132 KMeasured spectrum of CrA-9b Figure G1.
BT-SETTL models retrieved by specfit when using a gaussian prior of 𝜇 = . 𝜎 = . ( 𝑔 ) . This leads to only a slightly poorer fitthan log ( 𝑔 ) ∼ .
6, obtained with uniform priors. Favoured parameters are 𝑇 𝑒 = ( 𝑔 ) = . 𝑅 𝑏 = . 𝑅 𝐽 and 𝐴 𝑉 = . μ m)0.00.51.01.52.02.53.03.5 λ F λ ( W m − μ m − ) T e = 1655 KMeasured spectrum of CrA-9b0.8 1.3 1.8 2.3 2.8 3.3 3.8Wavelength ( μ m)0.00.51.01.52.02.5 λ F λ ( W m − μ m − ) T e = 1675 KMeasured spectrum of CrA-9b Figure G2.
BT-SETTL (top) and DRIFT-PHOENIX (bottom) models retrieved by specfit when fixing the photometric radius to 1.8 𝑅 𝐽 . Poor fits areobtained when fixing the photometric radius to the expected physical radius of a young Jovian planet with the absolute near-IR magnitudes of the companion.MNRAS000