Detectability of embedded protoplanets from hydrodynamical simulations
MMNRAS , 1–18 (2020) Preprint 13 January 2020 Compiled using MNRAS L A TEX style file v3.0
Detectability of embedded protoplanets fromhydrodynamical simulations
E. Sanchis, , (cid:63) G. Picogna, B. Ercolano , L. Testi , , and G. Rosotti , European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei M¨unchen, Germany Universit¨ats-Sternwarte, Ludwig-Maximilians-Universit¨at M¨unchen, Scheinerstrasse 1, D-81679 M¨unchen, Germany Excellence Cluster Origins, Boltzmannstrasse 2, D-85748 Garching bei M¨unchen, Germany INAF/Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, UK Leiden Observatory, Leiden University, P.O. Box , NL-
RA Leiden, the Netherlands
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We predict magnitudes for young planets embedded in transition discs, still affectedby extinction due to material in the disc. We focus on Jupiter-size planets at a latestage of their formation, when the planet has carved a deep gap in the gas and dustdistributions and the disc starts being transparent to the planet flux in the infrared(IR). Column densities are estimated by means of three-dimensional hydrodynamicalmodels, performed for several planet masses. Expected magnitudes are obtained byusing typical extinction properties of the disc material and evolutionary models of giantplanets. For the simulated cases located at . in a disc with local unperturbedsurface density of
127 g · cm − , a M J planet is highly extincted in J -, H - and K -bands, with predicted absolute magnitudes ≥
50 mag . In L - and M -bands extinctiondecreases, with planet magnitudes between and
35 mag . In the N -band, due to thesilicate feature on the dust opacities, the expected magnitude increases to ∼
40 mag .For a M J planet, the magnitudes in J -, H - and K -bands are above
22 mag , while for L -, M - and N -bands the planet magnitudes are between and
20 mag . For the M J planet, extinction does not play a role in any IR band, due to its ability to open deepgaps. Contrast curves are derived for the transition discs in CQ Tau, PDS , HL Tau,TW Hya and HD . Planet mass upper-limits are estimated for the known gapsin the last two systems. Key words:
Protoplanetary discs – Hydrodynamics – planet–disc interactions –planets and satellites: detection – infrared: planetary systems
For the last two decades a large number of exoplanet de-tections have remarkably expanded and shaped the pre-vailing planet formation theories. Most of these discoverieshave been accomplished using indirect techniques, althoughin a few cases detections were achieved using direct imag-ing (Quanz et al. 2015; Reggiani et al. 2018; Keppler et al.2018). Direct observations are only possible for systems closeenough to us, with planets far from their host stars (Rameauet al. 2015). When searching for young planets of only afew
Myr an additional problem arises: while a young planetmight itself be bright enough to be detected, it will still beembedded in the primordial disc, and thus hidden behindlarge columns of dust and gas. (cid:63)
E-mail: [email protected]
Studying the interaction between the disc material andthe young planet during its formation is crucial to under-stand the ongoing processes of planet formation and its fur-ther evolution. After the initial stages of planet formation,the protoplanet is likely accreting material and still sur-rounded by gas and dust. Several processes –mainly internaland/or external photo-evaporation, accretion onto the cen-tral star and planet, magneto-hydrodynamical winds (see re-view by Ercolano & Pascucci 2017) continuously reduce thedisc material until its complete dispersal. The typical life-times for discs can vary considerably and are uncertain, ingeneral the inner disc (within a fraction of AU from the cen-tral star) is expected to disperse within a few
Myr (Hern´an-dez et al. 2008), even though there is potential evidence forreplenishment of the inner disc over longer timescales (e.g.Beccari et al. 2010; Scicluna et al. 2014). The dissipationtimescales of the outer disc, which is more relevant for the © a r X i v : . [ a s t r o - ph . E P ] J a n Sanchis et al. direct detectability of young protoplanets, is much more un-certain. Initial Atacama Large Millimeter/submillimeter Ar-ray (ALMA) surveys suggest that the depletion of the outerdisc may also proceed on a similar timescale, following a sim-ple estimate of disc masses based on (sub)-millimetre con-tinuum emission from large dust grains (Ansdell et al. 2017).More detailed studies seem to imply that gas and dust deple-tion may be substantial even at young ages (Miotello et al.2017; Manara et al. 2018). As the disc keeps losing material,the extinction is reduced proportionally. If dust extinctionhas decreased enough, a detection of an embedded planetmay be possible with state-of-the-art facilities observing atIR wavelengths, where the planet spectrum peaks.Consequently, detections of protoplanets embedded indiscs depend on the properties of the planet, its immedi-ate surroundings, and also on the upper atmospheric lay-ers of the disc. The search for indirect detections in (sub)-millimetre observations has been pursued during the pastyears (e.g., Pi´etu et al. 2006; Brown et al. 2008), more in-tensively once ALMA started operating (ALMA Partner-ship et al. 2015). Substructure like cavities, gaps and spi-rals could be first observed at these wavelengths, suggestingplanet-disc interaction as a plausible cause of such features(e.g. Paardekooper & Mellema 2006; Rice et al. 2006). TheDSHARP large program (Andrews et al. 2018) has con-firmed that substructure is ubiquitous in large discs whenobserved with enough resolution, although these disc fea-tures do not necessarily confirm the presence of planets. Thefirst indirect detection of planets at these wavelengths wasachieved from detailed analyses of the gas kinematics in theHD disc (Teague et al. 2018; Pinte et al. 2018).While these indirect detections with ALMA and other(sub)-millimetre facilities tantalise evidence for young, andin some cases, massive planets in discs, the interpretation isnot unique. Direct detection of young planet candidates is re-quired to confirm their presence in discs. Additionally, directdetection from IR and spectroscopy is crucial to further char-acterise the planet properties (e.g. its atmosphere). Severalattempts have been made to detect directly young planetcandidates in discs, by means of near and mid-IR high con-trast imaging, but the vast majority of these efforts resultedin no detections. Testi et al. (2015) searched for planets inHL Tau in L -band, without any point-sources detected butsetting upper limits of planet masses at the rings location.Much more stringent upper limits were set for the TW Hyasystem by Ruane et al. (2017) with the Keck/NIRC2 in-strument, and for the HD system by Guidi et al.(2018). In only a few occasions, point-like sources have beenclaimed as detections in other protoplanetary discs. Reggianiet al. (2014) and Quanz et al. (2015) announced direct evi-dence (as a point-sources) of protoplanets embedded in theHD 169142 and HD 100546 at L - and M (cid:48) -bands respectively,using the NaCo instrument at the Paranal Observatory ofthe European Southern Observatory. While apparently con-vincing, both detections have been disputed and further con-firmation is still pending. More recently, an additional candi-date, still requiring confirmation, has been detected by Reg-giani et al. (2018) in the spiral arm of MWC 758 in L -band.To date, the most convincing direct detection of a youngplanet in a protoplanetary disc, is the multi-wavelength de-tection of PDS b, confirmed using multiple epoch data (Keppler et al. 2018); the system also includes an additionalplanet, PDS c (Haffert et al. 2019).Hydrodynamical (HD) simulations of planet-disc inter-actions have greatly improved our theoretical understandingof these systems. Kley (1999) and Bryden et al. (1999) sim-ulated a planet embedded in a disc in D, showing the for-mation of gaps from the planet-disc interaction and derivingcrucial properties of protoplanetary discs as mass accretionand viscosity. Many other studies followed, focusing on thecharacterisation of geometrical properties and the deriva-tion of important disc evolution parameters (like the viscos-ity, e.g. Szul´agyi et al. 2014). The implementation of nestedgrids alleviated the resolution limitation and computationtimes, allowing for the first D simulations to be performed(see D’Angelo et al. 2003).The detectability of a protoplanet at different wave-lengths can also been inferred from HD-simulations. Indeed,Wolf & D’Angelo (2005) performed mock observations to in-fer planet signatures that could be detected using ALMA,while Zhu (2015) focused on the detectability in IR-bandsof the circumplanetary disc (CPD) around highly accret-ing planets. Other indirect observable signatures due to thepresence of planets have been also investigated: gap openingeffects at various wavelengths (e.g. Dong et al. 2015a; Rosottiet al. 2016; Dipierro & Laibe 2017; Dong & Fung 2017;Jang-Condell 2017), disc inclination effects (Jang-Condell& Turner 2013), spiral arms in scattered light (Dong et al.2015b; Fung & Dong 2015; Juh´asz et al. 2015; Juh´asz &Rosotti 2018), planet shadowing (Jang-Condell 2009), andeven the effect of migrating planets (Meru et al. 2019; Nazariet al. 2019).This study focuses on the early evolutionary stages,when the accretion of disc material onto the planet mightstill have an incidence on the total planet brightness. Weperformed 3D HD simulations, using high resolution nestedgrids in the planet surroundings. Column densities and ex-tinction coefficients at the planet location are derived fromthe simulations, in order to infer planet magnitudes in J -, H -, K -, L -, M - and N -bands. This model can be used toguide future direct imaging observations of young planetsembedded in protoplanetary discs.The work is organised as follows: in Section 2 we de-scribe the simulations set-up and the model. The results ofthe different simulated systems are presented in Section 3.The application of the model to known protoplanetary discsis discussed in Section 4, and the main implications of thiswork are summarised in Section 5. The HD simulations were performed with the PLUTO code(Mignone et al. 2010, 2012), a 3-dimensional grid-code de-signed for astrophysical fluid dynamics. The details of thesimulations set-up, the description of the planet flux modelused (2.2) and the derivation of the expected planet magni-tudes (2.3) are presented in this section.
We modelled a gaseous disc in hydrostatic equilibriumaround a central star with a protoplanet on a fixed orbit.
MNRAS000
MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations The disc evolution is determined by the α -viscosity prescrip-tion (Shakura & Sunyaev 1973). The disc is considered to belocally isothermal, for which the Equation-of-State (EoS) isdescribed as: P = n · k B · T = ρ m u · µ · k B · T (1)where P is the pressure, n the total particle number density, k B the Boltzmann constant, T the temperature, ρ the gasdensity, m u the atomic mass unit and µ the mean molecularweight. The temperature in the disc varies only radially: T ( R ) = T (cid:18) RR (cid:19) q (2)with R being the radius in cylindrical coordinates, R = . AU, T = K and q the temperature exponent factor, setto − . The adopted density distribution of the protoplane-tary disc (as in Nelson et al. 2013) is described as: ρ ( r , θ ) = ρ (cid:18) RR (cid:19) p exp (cid:20) GM (cid:63) c · (cid:16) r − R (cid:17)(cid:21) (3)where r refers to the radial distance from the centre in spher-ical coordinates, θ the polar angle (thus R = r · cos ( θ ) ), ρ the density at the planet location, p is the density expo-nent factor, in our case with value − . , G is the universalgravitational constant, M (cid:63) the mass of the central star, and c iso the isothermal speed of sound. From this equation, thesurface density Σ scales radially as: Σ ( R ) = Σ · (cid:18) RR (cid:19) − p (cid:48) (4)with a power law with index p (cid:48) = − , . The planet is in-cluded as a modification in the gravitational potential ofthe central star in the vicinity of the planet location, whichis kept fixed. The gravitational potential φ considered in thesimulations is: φ = φ (cid:63) + φ pl + φ ind (5) φ (cid:63) is the term due to the star, φ pl the planet potential and φ ind accounts for the effect of the planet potential onto thecentral star. In the cells closest to the planet (cells at adistance to the planet lower that the smoothing length d rsm ), φ pl is introduced with a cubic expansion as in Klahr & Kley(2006) to avoid singularities at the planet location: φ pl ( d < d rsm ) = − GM pl d (cid:20)(cid:16) dd rsm (cid:17) − (cid:16) dd rsm (cid:17) + (cid:16) dd rsm (cid:17)(cid:21) (6)with d referring to the distance between cell and planet.We set d rsm to be . R Hill in the first 4 simulations. Intwo additional runs we decreased this value and doubled thegrid resolution in order to investigate the behaviour of thedensity in the cells close to the planet surface. In Table 1 allthe values of d rsm used in the simulations are shown.The gas rotational speed ( Ω ) is sub-keplerian; the addi-tional terms arise from the force balance equations in radialand vertical directions (see e.g. Nelson et al. 2013). It isdescribed by: Ω ( R , z ) = Ω K (cid:20) ( p + q ) (cid:16) hR (cid:17) + ( + q ) − qR (cid:112) R + z (cid:21) / (7)with z being the vertical coordinate, Ω k the keplerian orbitalspeed, and h the vertical scale-height of the disc. The simulated protoplanetary discs had a stellar massof . M (cid:12) and a surface density at . AU of 127 g · cm − ,consistent with estimates of the Minimum Mass Solar Neb-ula (Hayashi 1981) and the densest discs in the star formingregions in the Solar neighbourhood. (Andrews 2015; Tazzariet al. 2017). The model can be adapted to systems withdifferent characteristics by re-scaling our results to differentsurface densities and other parameters, as shown in section4. We modelled three different planetary masses ( , and M J ) embedded in a viscous disc with α = . . A fourthinviscid run with a M J was performed in order to study theeffect of viscosity on planet detectability, since recent studieswith non-ideal magneto-hydrodynamical effects have shownthat the disc can be laminar at the mid-plane (Paardekooper2017). Besides, two additional runs with doubled grid resolu-tion were done for the and M J cases, in order to improveour understanding of the disc-planet interaction at the up-per atmospheric layers of the planet. In total 6 simulationswere carried out, summarised in Table 1. During the first orbits the planetary mass was raised until its final value ( , or M J ) following a sinusoidal function, in order to preventstrong disturbances in the disc. The simulations ran for orbits, once a steady-state of the system is reached.The resolution at the vicinity of the planet is crucial forour study in order to obtain realistic infrared (IR) opticaldepths due to the disc material close to the planet surface.To fulfil this, we used a 3-level nested grid with a maximumresolution of . R Hill . For the , and M J planets, andconsidering planetary radii of . , . and . R J fromthe 1 Myr old hot-start models of Spiegel & Burrows (2012)(details of these models in Section 2.2), the maximum resolu-tion corresponds to . , . and . planetary radii respec-tively. In the 2 additional runs with doubled resolution overthe entire grid, the smallest cell-size was set to . R Hill .To test the accuracy of our simulations, we inspected thegas streamlines in the vicinity of the planet (Section 3.1.3).The grid is centred at the star location, thus close to theplanet the grid appears to be Cartesian. As suggested byOrmel et al. (2015a,b), the grid describes the CPD correctlyif the gas streamlines form enclosed circular orbits aroundthe planet location, which is confirmed in our simulations(Figure 1).To save computational time, we assumed the disc tobe symmetric with respect to the mid-plane. The simulatedrange for θ goes from the mid-plane up to ◦ , adequate for aproper representation of the disc dynamics and the columndensity derivation. At the disc upper layers, density has de-creased ≥ - orders of magnitude, and the contribution tothe column densities of these upper layers is negligible. Thedetails for the 3 grid regions are summarised in Table 2. The total planet emission is considered as a combination ofits intrinsic and accretion flux components. At the wave-lengths studied, the intrinsic flux of the planet is expectedto dominate, except in those cases with very high accretionrates onto the planet.The intrinsic component of the planet flux is derivedfrom the evolutionary models of Spiegel & Burrows (2012).These models provide the absolute magnitudes in J -, H -, K -, MNRAS , 1–18 (2020)
Sanchis et al.
Table 1.
Set-up parameters of the simulations and inferred column mass densities σ . The disc aspect ratio H , defined as H = h / R wasset to . in all the simulations. The columns refer to: the α -viscosity parameter for a viscously evolving disc; cell size at planet vicinity,given in Hill radii ( R Hill ) and in planet radii ( R pl ); smoothing and accretion radii; and predicted column mass densities in g · cm − . σ isobtained integrating over every cell above the planet except the cells within the d rsm . The uncertainty is computed from the dispersionof the σ value in the last orbits. Planet radii are taken from the evolutionary models of giant planets of Spiegel & Burrows (2012).Run α Cell size [ R Hill ] Cell size [ R pl ] d rsm [ R Hill ] r sink [ R Hill ] σ M J .
003 0 .
02 7 . .
10 0 .
07 2 . ± . M J .
003 0 .
02 9 . .
10 0 .
07 0 . ± . M J .
003 0 .
02 11 . .
10 0 .
07 0 . ± . M J , inviscid disc ∼ .
02 7 . .
10 0 .
07 1 . ± . M J , doubled resol. .
003 0 .
01 3 . .
03 0 .
03 2 . ± . M J , doubled resol. .
003 0 .
01 5 . . - .
05 0 .
03 0 . ± . Table 2.
Resolution and extension (as of cells) for each coordinate (radial distance R , polar θ and azimuthal φ angles) of our 3-levelsgrid. There is no low-resolution level for the θ coordinate.Level 1 ( hi-res ) Level 2 ( mid-res ) Level 3 ( low-res )Coord. Resolution . R Hill
128 0 . R Hill
64 128 θ . R Hill
64 0 . R Hill - φ . R Hill
128 0 . R Hill
128 128
Figure 1.
Density map at the mid-plane near a M J planet. Thegas streamlines are plotted on top, showing the circular motionof the gas around the planet. A value of 1 in radial code units isequivalent to . . The colour scale shows the logarithm of thedensity, in g · cm − . L -, M -, N - IR bands for a range of planet masses and as afunction of age, up to M y r . Following their nomenclaturewe refer to hot-start and cold-start models to the cases of aplanet fully formed via disc instability or via core accretionrespectively. In a more realistic scenario a planet would becharacterised by an intermediate solution.The total accretion luminosity is given by (Frank et al.1985): L acc (cid:39) GM pl (cid:219) M acc R acc (8)where M pl is the planetary mass, (cid:219) M acc is the mass accre-tion rate onto the planet and R acc the accretion radius, thedistance to the planet at which accretion shocks occur (Hart-mann 1998). R acc is typically - times the planetary radius, we have chosen R acc ≡ R pl for the estimation of the accretionluminosity. The accretion shocks are not resolved in the sim-ulations, this is out of the scope of this work. Nevertheless,the computation of (cid:219) M acc from the simulations (explained inlast paragraph of this section) are independent of R acc . Theaccretion flux considered in this model accounts only for theaccretion shocks’ irradiation. Flux irradiated by the CPD isnot taken into account.The continuum emission of the accretion shocks is attemperatures ( ≡ T acc ) of the order of ∼ K (Hartmannet al. 2016). To obtain the accretion flux in each band, we ap-proximate the shocks’ emission to a black body (Mendigut´ıaet al. 2011) that emits at a radius R acc . The surface areacovered by the shocks is a fraction (defined as b acc ) of thespherical surface with same radius. In this model, the accre-tion flux in a given band ( F bandacc ) is computed as a fraction b acc of the flux within the same band of a spherical blackbody ( F bandbb ) with temperature T acc and radius R acc : F bandacc = b acc · F bandbb (9)The factor b acc can be estimated as the fraction be-tween the total accretion luminosity L acc (from Eq. 8) andthe bolometric luminosity of the same black body: b acc = L totacc L totbb (10)We tested various values of T acc , the results shown along thiswork are obtained using a shock temperature of . Asdiscussed for the band fluxes results (Section 3.2.1), assum-ing lower T acc does not vary the fluxes results significantly.The gas accretion onto the planet is modelled follow-ing the prescription in Kley (1999), also D¨urmann & Kley(2015). At each time-step, a fraction f acc · ∆ t · Ω of the gas is re-moved from the cells enclosed by a sphere of radius r sink cen-tred at the planet (values for r sink in Table 1). This methodmimics the direct accretion onto the planet surface. The val-ues of r sink have been chosen to guarantee convergence; asshown in Tanigawa & Watanabe (2002), the estimated ac-cretion rates converge to a stable value if r sink (cid:46) . R Hill .This is independent of the f acc value, that is set to . The MNRAS000
Density map at the mid-plane near a M J planet. Thegas streamlines are plotted on top, showing the circular motionof the gas around the planet. A value of 1 in radial code units isequivalent to . . The colour scale shows the logarithm of thedensity, in g · cm − . L -, M -, N - IR bands for a range of planet masses and as afunction of age, up to M y r . Following their nomenclaturewe refer to hot-start and cold-start models to the cases of aplanet fully formed via disc instability or via core accretionrespectively. In a more realistic scenario a planet would becharacterised by an intermediate solution.The total accretion luminosity is given by (Frank et al.1985): L acc (cid:39) GM pl (cid:219) M acc R acc (8)where M pl is the planetary mass, (cid:219) M acc is the mass accre-tion rate onto the planet and R acc the accretion radius, thedistance to the planet at which accretion shocks occur (Hart-mann 1998). R acc is typically - times the planetary radius, we have chosen R acc ≡ R pl for the estimation of the accretionluminosity. The accretion shocks are not resolved in the sim-ulations, this is out of the scope of this work. Nevertheless,the computation of (cid:219) M acc from the simulations (explained inlast paragraph of this section) are independent of R acc . Theaccretion flux considered in this model accounts only for theaccretion shocks’ irradiation. Flux irradiated by the CPD isnot taken into account.The continuum emission of the accretion shocks is attemperatures ( ≡ T acc ) of the order of ∼ K (Hartmannet al. 2016). To obtain the accretion flux in each band, we ap-proximate the shocks’ emission to a black body (Mendigut´ıaet al. 2011) that emits at a radius R acc . The surface areacovered by the shocks is a fraction (defined as b acc ) of thespherical surface with same radius. In this model, the accre-tion flux in a given band ( F bandacc ) is computed as a fraction b acc of the flux within the same band of a spherical blackbody ( F bandbb ) with temperature T acc and radius R acc : F bandacc = b acc · F bandbb (9)The factor b acc can be estimated as the fraction be-tween the total accretion luminosity L acc (from Eq. 8) andthe bolometric luminosity of the same black body: b acc = L totacc L totbb (10)We tested various values of T acc , the results shown along thiswork are obtained using a shock temperature of . Asdiscussed for the band fluxes results (Section 3.2.1), assum-ing lower T acc does not vary the fluxes results significantly.The gas accretion onto the planet is modelled follow-ing the prescription in Kley (1999), also D¨urmann & Kley(2015). At each time-step, a fraction f acc · ∆ t · Ω of the gas is re-moved from the cells enclosed by a sphere of radius r sink cen-tred at the planet (values for r sink in Table 1). This methodmimics the direct accretion onto the planet surface. The val-ues of r sink have been chosen to guarantee convergence; asshown in Tanigawa & Watanabe (2002), the estimated ac-cretion rates converge to a stable value if r sink (cid:46) . R Hill .This is independent of the f acc value, that is set to . The MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations mass accreted is removed from the computational domaininstead of being added to the planet mass (accreted massis negligible compared to the total planet mass, about − times lower). The extinction in our model is due to the disc materialaround the planet. We assume that there are no additionalastronomical objects of significant brightness or size betweenthe studied planet and the observer. Extinction in the V -band is derived from the column mass density obtained fromthe simulations. Using the magnitude-flux conversion for-mula from (G¨uver & ¨Ozel 2009), assuming constant gas-to-dust ratio of and a molecular weight of . : A V [ mag ] = N H [ atoms · cm − ] . · [ atoms · cm − · mag − ] (11)The above relation was inferred from observations andit applies to an averaged interstellar medium (ISM) in theMilky Way. While the gas-to-dust ratio used is a good firstapproximation, variations may be expected, particularly inthe vicinity of the planet. Recent surveys in close star-forming regions indicate that this ratio might indeed be dif-ferent in many discs (Miotello et al. 2017).From the inferred A V , the extinction coefficients ( A band )and optical depths ( τ band ) are obtained using the diffuse ISMextinction curves of Cardelli et al. (1989) for J -, H -, K -bands,and Chiar & Tielens (2006) (which accounts for the silicatefeature around µ m ) for L -, M - and N -bands. The expectedfluxes are then given by: F bandexpected = F bandpl · e − τ band (12)with F bandpl being the total planet flux in that band. Theresulting F bandexpected is the expected band flux that we wouldobserve for a planet embedded in a protoplanetary disc withongoing accretion onto the planet. In this section we first present the results from the HD sim-ulations (subsection 3.1), followed by the results from ourmodel (subsection 3.2)
Jupiter-size planets embedded in a disc generally carve a gapafter several orbits (Bryden et al. 1999). In our simulations,this can be seen in 2D density maps of the disc as a functionof time (Figure 2, for the M J case). For more massive plan-ets, the gap carving process is faster, as one would expectfrom planet formation theory.The gap opened by each planet can be compared to discmodels to verify the quality of our simulations. We tested ourresulting surface density profiles with an analytical model forgaps in protoplanetary discs, as described in Duffell (2015).An algebraic solution of the gap profiles is presented in thatwork, together with the derivation of a formula for the gap depth. In Figure 3 we show the azimuthally averaged sur-face density radial profiles relative to the unperturbed sur-face density ( Σ ) for the , and M J simulated planetsin a viscous disc. The solid lines represent the profiles aftera steady-state is reached, and the dashed lines denote thesurface density after the first orbits. The predicted gapdepths from the model (Duffell 2015, equation 9) are shownin the figure as horizontal dash-dot lines.Our results for and M J planets are in very goodagreement with the analytical model. The gap for a M J planet is relatively deeper than the prediction from themodel. Nevertheless, the model by Duffell (2015) fails atreproducing the gap profile produced by planets with veryhigh masses, as discussed in that work. Therefore, we cantrust the quality of our simulations from the concordance atlow planet masses with analytical models. The column density is obtained by integrating the disc den-sity over the line-of-sight towards the planet. We consideredour system to be face-on, since is the most likely geometryfor a direct protoplanet detection. The density above theplanet for every simulation is shown in Figure 4, togetherwith the unperturbed (initial) density. The highest densitiescorrespond to the M J simulation. For more massive plan-ets, the disc densities are lower, since the planet carves adeeper gap. In the inviscid case, the carved gap can not berefilled with adjacent material, resulting in lower densitiescompared to the viscously evolving case.The planet and its atmosphere are not resolved in thesesimulations, therefore the first cell is assumed to be theplanet outer radius. The sharp peak from the vertical den-sity profiles extends over the first 2-3 cells. This is expectedto be a combination of several effects, mainly an artefactof the simulations due to the potential smoothing (Eq. 6),which affects every cell within d rsm . However we cannot ex-clude the possibility that a fraction of the peak might alsobe the real density stratification of the material at the layersclosest to the planet. Additionally, this over-density is alsoaltered by the accretion radius (see 2.2), and limited by thegrid resolution. A fully resolved CPD would be necessary todisentangle between the different causes. However, we testedwhether the extension of the peak is an artefact due to thepotential softening and the mentioned resolution limitations.To this aim we performed 2 additional simulations with dou-bled resolution over the entire grid. We also decreased bothaccretion and smoothing radii to . - . R Hill . The val-ues of the main parameters for the doubled resolution runsare summarised in Table 1. This test was done for the discswith and M J planets. The grid from the last snapshot ofthe original simulations were readjusted to the new resolu-tion, and we let the system evolve until a new steady-statewas reached ( ≤ orbits needed in both cases). The verticaldensity at the closest cells above the planet for both origi-nal and doubled resolution for the M J case are shown inFigure 5. The vertical grid-lines represent . R Hill , and thecoloured lines illustrate the d rsm in both original and dou-bled resolution runs. The figure shows that the over-densitywith doubled resolution spans approximately half the orig-inal case. This is in accordance with what we expect if theover-density is due to the smoothing within d rsm . If we were MNRAS , 1–18 (2020)
Sanchis et al.
Figure 2.
Evolution in time of the 2D density map at the mid-plane for the viscous disc with an embedded M J planet. A value of in radial code units is equivalent to . (Jupiter semi-major axis). Density is represented in logarithmic scale, with values in g · cm − .From top left to bottom right, each snapshot represents the density ρ of the disc after , , , , and orbits. Figure 3.
Surface density radial profile for , and M J planetsafter the simulations reach a steady-state, shown as solid lines.The dash-dot lines are the respective predicted gap-depths, de-rived from an analytical model for gaps in protoplanetary discs(Duffell 2015). The dashed lines represent the surface density af-ter orbits. A value of in radial code units is equivalent to . . able to completely remove the potential smoothing we wouldonly see a peak inside the planet radius, i.e. within the firstcell. From the results of this test, we consider the columnmass density σ as the integrated density for all the cellsabove the smoothing radius. The resulting σ for each of thesimulated systems are included in Table 1. The uncertaintyconsidered is the dispersion of the column mass density forthe last orbits of each simulation. The predicted magni-tudes for and M J planets are derived using the σ val-ues from the doubled resolution runs, since in these cases Figure 4.
Vertical density above the planet for each simula-tion. Horizontal axis represents the height from the mid-planein Jupiter semi-major axis ( a J ); vertical axis shows density inlogarithmic scale, in g · cm − units. the planet-disc interaction is represented more accurately.The values for the original and doubled resolution cases arewithin their respective uncertainty: for a M J we obtained . ± . and . ± . · cm − , while for the M J runs the col-umn mass densities were . ± . and . ± .
01 g · cm − respectively. The prescription used to derive the mass accretion ratesfrom the simulations was described in Subsection 2.2. Thefinal value of (cid:219) M acc for each simulated planet is the opti-mal value of parameter c when fitting the evolution of (cid:219) M acc MNRAS000
01 g · cm − respectively. The prescription used to derive the mass accretion ratesfrom the simulations was described in Subsection 2.2. Thefinal value of (cid:219) M acc for each simulated planet is the opti-mal value of parameter c when fitting the evolution of (cid:219) M acc MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations Figure 5.
Vertical density profile close to the M J planet, com-paring the doubled resolution (blue line) to the original case (red).The vertical grid spacing is . R Hill , equivalent to the cell-sizefor the doubled resolution run, and half the cell-size of the origi-nal run. The red and blue dashed lines represent the smoothingradii for the original and the new run respectively. in the last orbits by an exponential function defined as f ( x ) = a · e − b · x + c . The resulting accretion rates are of theorder of − M (cid:12) · yr − , summarised in Table 3. The high-est (cid:219) M acc is obtained for a system with M J planet. Thisis somewhat counter-intuitive, as one may expect the leastmassive planet to have the lowest accretion rate, and themost massive planet to have the highest. There are two ef-fects to consider, the viscosity that refills the gap and thegravitational force of the planet: a more massive planet cre-ates a stronger gravitational potential, a larger CPD andhas a larger accretion rate, on the other hand if the discevolves viscously it can replenish the accreted material withnew material, thus keeping a higher (cid:219) M acc value. Beyond acertain planet mass, the gravitational potential clears largeregions quicker, which limits the refill of the gap materialthus ultimately reducing the (cid:219) M acc . Our results indicate thatthe mass of the planet for which this occurs is between and M J . The different accretion rates obtained for the M J planet in a viscous and an inviscid disc can be under-stood by considering that an inviscid disc cannot replenishthe gap opened by a planet, and consequently there is lessmaterial to feed the CPD.Care should be taken when considering these resultsbecause of the limitations posed by isothermal simulations,which do not account for accretion heating in the vicinityof the planet. This results in an overestimation of the ac-cretion rate. These simulations provide an upper limit ofthe expected accretion rates, and this should be taken intoaccount when the model is applied to real systems. Gas streamlines provide useful insights on the disc behaviourand the planet-disc interaction. In the close-up top view ofthe system (Figure 1), the gas streamlines at the system mid-plane are plotted as vectors. As expected from the accretionmodels, a CPD is formed and the gas motion is concentric tothe planet. The gas streamlines follow closed trajectories inthe regions nearest to the planet, which indicates that the
Figure 6.
Edge on view density map with gas streamlines aroundthe planet. An important part of the gas is being accreted ver-tically. A value of 1 in radial code units is equivalent to . .The scale represents the logarithm of the density, in g · cm − . Table 3.
Mass accretion rates and bolometric fluxes for everysimulated system. Fluxes are expressed in [ W · m − ]. M pl (cid:219) M acc [ M (cid:12) · yr − ] F Hotintr F Hotacc F Coldintr F Coldacc M J . e − . e . e . e . e M J . e − . e . e . e . e M J . e − . e . e . e . e M J inviscid . e − . e . e . e . e grid does describe the accreting planet with a CPD accu-rately (as discussed in Ormel et al. 2015a).Accretion onto the planet occurs not only in the orbitalplane from the CPD but also vertically (see e.g. Szul´agyiet al. 2016). In our simulated systems, the gas is indeedfalling onto the planet from its pole or with small inclinationangle from the vertical. Figure 6 shows that a high fractionof the material is being accreted vertically onto the M J planet. Accretion and intrinsic bolometric fluxes for each planet con-sidered in this work are shown in Table 3. Mass accretionrates are also included in the table. The difference in F acc between hot and cold-start models arises from the differentplanet radius R pl of each model (taken from Spiegel & Bur-rows 2012), since the R acc used to compute the accretionflux (Equation 8) is assumed to be ≡ R pl . The bolometricaccretion flux is higher than the planet’s intrinsic flux for allcases. Nevertheless, the accretion flux (which peaks at . µ m , with T eff ∼ ) is in most cases lower in the IRbands considered than the intrinsic planet flux, whose spec-trum peaks in the IR ( ∼ - µ m ). This can be seen in theleft panel of Figure 7. The figure shows the accretion andintrinsic fluxes of planets with , and M J for hot and cold-start planet models. Intrinsic fluxes are considerably lowerfor cold-start planets than for the hot models. The accretioncontribution can be significant, especially in the cold-start cases in J -, H -, K - and L -bands. Reducing the T acc shiftsthe accretion maximum to longer wavelengths, but it highlydecreases its bolometric value, and consequently the overallpicture does not vary significantly. MNRAS , 1–18 (2020)
Sanchis et al.
Figure 7.
Intrinsic ( F intr ) and accretion fluxes ( F acc ) for , and M J planets. The results for planets with a hot-start model are shownon the left panel, while cold-start models on the right. The intrinsic fluxes are plotted at the central wavelength of each band. These results indicate that radiation from accretionshocks near the planet is an important factor of the planetflux. Nevertheless, it is worth keeping in mind that these re-sults are for accretion rates inferred from isothermal simula-tions, which are generally overestimated. For more realisticaccretion rates at this stage ( ∼ − M (cid:12) · yr − ), intrinsic fluxdominates at these wavelengths. When scaling our models tolarger distances from the host star, accretion rates are highlyreduced due to the scaling for lower disc densities. Conse-quently intrinsic fluxes dominate in planets further out inthe disc for both hot and cold-start scenarios. From the column mass densities we derived extinction co-efficients for each planet in J - H -, K -, L -, M - and N -bands.The inferred coefficients are included in Table 4, and plot-ted on the top left panel in Figure 8. Extinction coefficientsat wavelengths (cid:46) µ m are extremely high for and M J .The effect of extinction decreases for longer wavelengths,but it raises up again at ∼ - µ m due to the silicate featurepresent in the diffuse ISM. For a M J , extinction coefficientsare very low in every IR band. This is due to its ability toopen a gap in the disc very effectively, and consequently, discdensity around the planet and the inferred column densityare very low.The derived magnitudes for each simulated planet as afunction of wavelength are shown on the top right panel inFigure 8. For each planet, the upper and lower curves rep-resent its hot and cold models. The results include extinc-tion from the disc material and radiation from the accretionshocks near the planet. The magnitudes for every IR-bandand planet at . are summarised in Table 4. For eachsimulated planet, the rows in the table show: the absolutemagnitude of the planet including the contribution from ac-cretion; the extinction coefficients due to the disc material,and the predicted absolute magnitude including extinctioneffects. The predicted magnitudes for a given planet decreaseas a function of wavelength, except at ∼ µ m due to thesilicate feature in dust grain opacities. The curves are moreflattened for more massive planets. This is due to the moreefficient depletion of material in the planet vicinity, whichyields lower column densities, hence lower extinction. Forthe most massive planet considered ( M J ), the gap clearingis extremely effective and extinction in any IR band is neg-ligible. In the inviscid scenario with a M J , the gap can notbe replenished efficiently compared to the viscously evolvingcase. This results on a gap with lower density and extinc-tion. At (cid:46) µ m , - M J planets are completely obscuredby the disc material (with extinction above 30 and 15 mag respectively). These results are obtained for an unperturbedsurface density of Σ =
127 g · cm − at . .We have used the ISM law to estimate the extinc-tion under the assumption that mostly small grains will bepresent in the disc atmosphere above the planet. The actualvalue of extinction depends on the assumption on the dustproperties. To investigate the implications of our assump-tion, we evaluated the impact of using different assump-tions on the dust composition and size distributions. Theresults for the M J viscously evolving case are shown inbottom panels of Figure 8. Two other dust models were in-vestigated: one model of grains with fractional abundancescomparable to the expected in protoplanetary discs mid-plane (Pollack et al. 1994) and grain population with num-ber density n ( a ) ∝ a − . (where a is the grain size) between . µ m < a < µ m (Tazzari et al. 2016, for details), and adust coagulation model for ice-coated silicate-graphite ag-gregates (type II grain mixing, see Ormel et al. 2009, 2011),applicable to dust in protoplanetary discs (extinction shownis for grain sizes a ∼ µ m ). The results of the ice silicategraphite model are in very good agreement with the ISMextinction used (in M -, and N -bands the diffuse ISM ex-tinction becomes larger). The other model provides similarresults in the J -, and H -bands, however, for longer wave-lengths extinctions are ∼ - times larger than for the diffuseISM. Thus, when considering dust with different properties MNRAS000
127 g · cm − at . .We have used the ISM law to estimate the extinc-tion under the assumption that mostly small grains will bepresent in the disc atmosphere above the planet. The actualvalue of extinction depends on the assumption on the dustproperties. To investigate the implications of our assump-tion, we evaluated the impact of using different assump-tions on the dust composition and size distributions. Theresults for the M J viscously evolving case are shown inbottom panels of Figure 8. Two other dust models were in-vestigated: one model of grains with fractional abundancescomparable to the expected in protoplanetary discs mid-plane (Pollack et al. 1994) and grain population with num-ber density n ( a ) ∝ a − . (where a is the grain size) between . µ m < a < µ m (Tazzari et al. 2016, for details), and adust coagulation model for ice-coated silicate-graphite ag-gregates (type II grain mixing, see Ormel et al. 2009, 2011),applicable to dust in protoplanetary discs (extinction shownis for grain sizes a ∼ µ m ). The results of the ice silicategraphite model are in very good agreement with the ISMextinction used (in M -, and N -bands the diffuse ISM ex-tinction becomes larger). The other model provides similarresults in the J -, and H -bands, however, for longer wave-lengths extinctions are ∼ - times larger than for the diffuseISM. Thus, when considering dust with different properties MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations Table 4.
Absolute magnitudes for planets at . to the host star, with masses: M J –viscous and inviscid scenarios– M J and M J . Mag pl is the total magnitude of the planet, including accretion flux; A band is the extinction coefficient in each band, Mag expected is themagnitude of the planet considering extinction due to disc material. All values are in mag . Hot-start planet
Cold-start planetJ H K L M N J H K L M N
Mag pl .
74 13 .
75 12 . .
72 13 .
81 13 .
84 13 .
56 12 .
65 11 . M J A band .
32 58 .
41 36 .
75 19 .
96 15 .
63 28 .
15 91 .
32 58 .
41 36 .
75 19 .
96 15 .
63 28 . expected .
07 72 .
16 49 .
67 32 .
36 26 .
96 38 .
16 105 .
04 72 .
21 50 .
59 33 .
52 28 .
28 40 . pl .
58 12 .
22 11 .
44 10 .
87 10 .
51 9 .
31 12 .
72 12 .
80 12 .
84 12 .
60 12 .
10 11 . M J A band .
33 16 .
84 10 .
60 5 .
84 4 .
57 8 .
24 26 .
33 16 .
84 10 .
60 5 .
84 4 .
57 8 . expected .
91 29 .
06 22 .
04 16 .
71 15 .
08 17 .
55 39 .
06 29 .
65 23 .
44 18 .
44 16 .
67 19 . pl .
09 10 .
21 9 .
48 9 .
00 9 .
29 8 .
33 12 .
61 12 .
69 12 .
73 11 .
88 11 .
62 11 . M J A band .
36 0 .
23 0 .
15 0 .
08 0 .
06 0 .
11 0 .
36 0 .
23 0 .
15 0 .
08 0 .
06 0 . expected .
45 10 .
45 9 .
63 9 .
08 9 .
35 8 .
45 12 .
97 12 .
92 12 .
88 11 .
96 11 .
68 11 . pl .
98 14 .
83 13 .
28 12 .
63 11 .
39 10 .
02 15 .
56 15 .
66 15 .
62 14 .
79 12 .
95 12 . M inviscid A band .
40 25 .
20 15 .
86 8 .
74 6 .
84 12 .
32 39 .
40 25 .
20 15 .
86 8 .
74 6 .
84 12 . expected .
38 40 .
03 29 .
14 21 .
37 18 .
24 22 .
35 54 .
96 40 .
86 31 .
48 23 .
52 19 .
79 24 . Figure 8.
Top panels: extinction coefficients with uncertainties (left), and predicted magnitudes for the simulated systems (right), bothas a function of wavelength. The predicted planet magnitudes are shown as an area delimited by the hot and cold planetary model.Bottom panels: extinction and predicted magnitudes of the M J viscous case using different dust grain models. The results for thevarious dust models are normalised at A V . For every panel, the vertical dotted lines represent (from left to right) the central wavelengthof J -, H -, K -, L -, M -, and N -bands.MNRAS , 1–18 (2020) Sanchis et al. (e.g. composition, size, level of processing), the resulting pre-dicted planet magnitudes might change due to the opacityvariations in IR wavelengths.On the other hand, the extinction is obtained assuminga gas-to-dust ratio of 100 along the disc. In the atmosphericlayers of the disc above the planet, this ratio might be largerdue to dust processing and settling. In this regard, our anal-ysis provides a conservative estimate of extinction, and thepresented results can be interpreted as the worst case sce-nario.
The simulations performed in this work are locally isother-mal, therefore the results can be scaled to account for differ-ent disc densities without altering the dynamics of the disc.We re-scaled column mass densities and accretion rates fordifferent planet positions and disc densities. The normalisa-tion of the surface density is readjusted in order to preservethe surface density profile Σ ( R ) as in Equation 4. In this waywe could extend the results of our model to study differentsystems. From dimensional analysis, the column mass den-sity σ is directly proportional to the surface density, whilethe accretion rate (cid:219) M acc is proportional to the density andinversely proportional to orbital time. Surface density Σ de-creases as R − . (Equation 4), while the density at the planetlocation ρ ∝ R − . (Equation 3). These relationships are usedto derive the scaling factors for σ and (cid:219) M acc for planets ata distance (cid:44) . . In case of re-normalising the surfacedensity (as done for real systems in Section 4) the ratio be-tween the unperturbed new and original surface densitiesat a fiducial distance multiplies the scaling factors of thecolumn mass density and the mass accretion rate.Scaling the distance to the central star would change thetemperature at the planet location. While this has no direct impact on the scaling applied to the (cid:219) M acc obtained from thesimulations (which has no explicit dependence on temper-ature), it would also change the disc aspect ratio H = h / R since discs are generally flared. This would have an indirect effect on (cid:219) M acc . In addition, the disc aspect ratio also deter-mines the gap opening planet mass and therefore the depthof the gap. We neglect these effects and remark that this isa limitation of our approach; the performed scaling providesa valuable understanding of how the planet location influ-ences accretion rates and densities, but it does not captureall possible effects.Tables with extinction coefficients and predicted mag-nitudes of the simulated systems at , , and forevery band are included in the Appendix A (Tables A1, A2,A3, A4). Figure 9 shows the expected absolute magnitudesof planets with (for both viscous and inviscid cases), and M J for different distances to the central star, in J -, H -, K -, L -, M , and N -bands. The coloured area of each planetrepresents the uncertainty associated to the column density.Extinction decreases for planets further out in the discdue to lower column densities. This behaviour is driven bythe surface density profile. Accretion flux is higher at shorterdistances but its effect is minor compared to extinction. Fur-ther out in the disc, the planet magnitude is dominated bythe intrinsic flux since the accretion drops with distance.The dispersion of the results is considerable, it can not beneglected as a source of indeterminacy in our results. Never- theless, large uncertainties are linked to very high values ofthe column density, and in such cases extinction is so strongthat the planet would be completely hidden.At shorter wavelengths ( J - H - and K -bands), the magni-tude of a M J planet is completely dominated by extinction:at the furthermost location considered,
100 AU , the extinc-tion is . , . , and . in each of these bands. For a M J planet the effect of extinction is lower but still large,e.g. at
100 AU the extinction in these bands is . , . , and . .The magnitude of a M J planet in the simulated discis barely affected by extinction in any band (the highest ex-tinction coefficient, A J , ranges between . to . between . and
100 AU ): since the planet is substantially more mas-sive, it has cleared almost all the material at the gap and itsvicinity, thus both column density and the extinction coeffi-cients are exceptionally low compared to the other simulatedplanets.In L -band, the disc material above the least massiveplanet causes an extinction of
10 mag at
20 AU , and is re-duced to . at
100 AU . A M J planet is less affectedby extinction, with A L of . and . at those distances.In the M -band, the extinction coefficients of all the planetsconsidered are the lowest out of all bands considered, nev-ertheless still considerable except for the M J planet. Forinstance, a M J planet at
100 AU would be extincted by . . In the N -band, extinction is increased due to thesilicate feature in the opacity of ISM dust: the extinctioncoefficients are almost 2 times larger than the coefficients inthe M -band.Our models can as well be used for different values ofthe stellar mass. Scaling our results for different M (cid:63) is espe-cially useful when applying our detectability model to realsystems, as discussed in the next section. The ratio M (cid:63) M pl can-not change in order to keep the dynamics of the systemvalid. The planet mass is re-scaled accordingly to keep thisratio constant. Since the planetary models from Spiegel &Burrows (2012) only provide data for planets with , , or M J masses, the planet intrinsic magnitudes were interpo-lated for the new M pl . Our results can be applied to real systems to study the de-tectability of embedded planets, proceeding as explained in3.2.3. In what follows we present the results of our modelfor the Class II discs of CQ Tau, PDS , HL Tau, TW Hyaand HD . For the last three systems, our results arecombined with contrast limits from previous IR observations(Testi et al. 2015; Ruane et al. 2017; Guidi et al. 2018). Theimprovement of this revision is the inclusion of the extinc-tion due to the disc material, and the emission from theshocks due to planet accretion.Additionally, to understand how likely would be to de-tect the simulated planets directly with ALMA, we esti-mated the M J planet and CPD fluxes at µ m wavelength.The expected planet flux is ∼ − mJy, below the ALMAsensitivity limit. For the CPD, simplified to a disc of R Hill radius centred at the planet and . R Hill high (i.e. the re-gion around the planet with a disc-shaped overdensity), weobtain a dust mass of M CPDdust ≈ . M ⊕ , which is compa- MNRAS000
100 AU would be extincted by . . In the N -band, extinction is increased due to thesilicate feature in the opacity of ISM dust: the extinctioncoefficients are almost 2 times larger than the coefficients inthe M -band.Our models can as well be used for different values ofthe stellar mass. Scaling our results for different M (cid:63) is espe-cially useful when applying our detectability model to realsystems, as discussed in the next section. The ratio M (cid:63) M pl can-not change in order to keep the dynamics of the systemvalid. The planet mass is re-scaled accordingly to keep thisratio constant. Since the planetary models from Spiegel &Burrows (2012) only provide data for planets with , , or M J masses, the planet intrinsic magnitudes were interpo-lated for the new M pl . Our results can be applied to real systems to study the de-tectability of embedded planets, proceeding as explained in3.2.3. In what follows we present the results of our modelfor the Class II discs of CQ Tau, PDS , HL Tau, TW Hyaand HD . For the last three systems, our results arecombined with contrast limits from previous IR observations(Testi et al. 2015; Ruane et al. 2017; Guidi et al. 2018). Theimprovement of this revision is the inclusion of the extinc-tion due to the disc material, and the emission from theshocks due to planet accretion.Additionally, to understand how likely would be to de-tect the simulated planets directly with ALMA, we esti-mated the M J planet and CPD fluxes at µ m wavelength.The expected planet flux is ∼ − mJy, below the ALMAsensitivity limit. For the CPD, simplified to a disc of R Hill radius centred at the planet and . R Hill high (i.e. the re-gion around the planet with a disc-shaped overdensity), weobtain a dust mass of M CPDdust ≈ . M ⊕ , which is compa- MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations Figure 9.
Absolute magnitudes at J -, H -, K -, L -, M , and N -bands for the simulated disc with an embedded planet ( , , M J , or a M J planet inviscid case) at different distances to the central star. The results are shown for a hot-start model of the formation scenario.The vertical axis for J and H -bands covers a wider range in order to include all the planets in the same panel. rable to the CPD measurements in PDS b (Isella et al.2019). Assuming a constant CPD temperature of
121 K , thecontinuum emission in ALMA Band 7 would be . mJyif emission is assumed optically thin, and . mJy if opti-cally thick. Thus, the CPD of the simulated disc could bedetected by ALMA observations with enough sensitivity. CQ Tau is a young star from the Taurus-Auriga region, spec-tral type A (Trotta et al. 2013) and M (cid:63) = . M (cid:12) (Testiet al. 2003). It has an estimated age of ∼ -
10 Myr (Chapillon et al. 2008), and very low disc mass, of the order of − - − M (cid:12) . A fiducial surface density of Σ = . · cm − at (a factor × . respect the simulated disc) was usedfor the re-normalisation to our unperturbed profile, derivedfrom ALMA observations (Ubeira Gabellini et al. 2019).We applied our models to investigate the effects of discextinction on potential planets embedded in the disc. Due tothe re-scaling with the stellar mass, the contrast curves arederived for planets with . , . and . M J . A distanceof . (Gaia Collaboration et al. 2018) was used forthe predicted contrast. Figure 10 shows the contrast of theplanets in the L -band as a function of distance to the central MNRAS , 1–18 (2020) Sanchis et al.
Figure 10.
Application of our model to planets with . , . and . M J masses embedded in the CQ Tau disc. The colouredlines represent the contrast in L -band of each planet as a functionof the distance to the central star placed at different distancesalong the disc. star; the planet contrast is shown relative to the stellar value.The coloured area for each planet represents the associateddispersion. The results in J -, H - and K -bands are includedin Appendix B.At a fixed distance, the more massive the planet, theless affected by extinction, since the planet is more effectiveat clearing the gap. Due to the very low surface density in-ferred from the ALMA observations, extinction in L -bandis only relevant for the lightest planets. At
20 AU , a . M J planet would have a contrast of .
16 mag and extinc-tion A L = .
34 mag ( A V = .
47 mag ). The contrast of a . M J planet is .
27 mag with an extinction of only . ( A V = .
60 mag ). The most massive planet ( . M J )is barely affected by extinction, its contrast is equivalent tothe case of a completely depleted disc, .
17 mag . PDS is a member of the Upper Centaurus-Lupus sub-group (at ∼
113 pc , Gaia Collaboration et al. 2018), with acentral star of . and mass . M (cid:12) (M¨uller et al. 2018).It is surrounded by a transition disc with estimated totaldisc mass of · − M (cid:12) . A first companion (PDS b) wasfound combining observations with VLT/SPHERE, VLT/-NaCo and Gemini/NICI at various epochs, detected as apoint-source in H -, K - and L -bands at a projected averagedseparation of . (Keppler et al. 2018). In J -band,PDS b could only be marginally detected when collapsingthe J - and H -band channels. Due to the high uncertainties, J -band magnitude was not given. Atmospheric modelling ofthe planet was used to constrain its properties (M¨uller et al.2018), with an estimated mass range from to M J .Recent H α line observations using VLT/MUSE con-firmed a σ detection from a second companion (PDS c)at
240 mas (Haffert et al. 2019). Dust continuum emission(likely from its CPD) has been also observed (Isella et al.2019). This second source is very close to an extended discfeature, consequently its photometry should be done withcaution. In Mesa et al. (2019), the planetary nature of thiscompanion has been confirmed, and absolute magnitudes in J -, H -, and K -bands could be inferred for two SPHEREepochs. The spectrum in the J -band is very faint, and indis-tinguishable from the adjacent disc feature, thus the J -bandmagnitude should be regarded as upper limit. The NaCo L -band map detected emission that is partly covered by thedisc, therefore its L -band magnitude should also be takenas an upper limit. Using various atmospheric models, Mesaet al. (2019) constrained the mass PDS c to be between . and . M J .Our models were re-scaled using a fiducial surface den-sity of Σ = . · cm − at (taking the unperturbedsurface density model with depletion factor δ disc = andgas-to-dust ratio of 100, from Keppler et al. 2018). This cor-responds to a surface density scale factor of × . with re-spect to the simulated disc. From the re-scaling, we obtainedcontrast curves of planets embedded in the PDS disc with . , . and . M J (Figure 11). The results show the ef-fect of a disc with very low surface density: extinction hasan incidence in J - and H -bands for . and . M J planetslocated within (cid:46)
40 AU . In the L -band extinction has only aminor effect on the lightest planet model at distances below
20 AU . From the assumed surface density profile, none of theplanetary companions would be affected by extinction dueto material from the protoplanetary disc in the IR bands.The observed contrast of the primary companion inthree bands is considerably higher than the value for themost massive planet of our models, thus setting a mass lowerlimit of . M J for PDS b. The second companion layson top of the . M J model in H -band, and above it in K -band. The redness of this source can explain the differencein the bands contrast. This reddening might be due to ma-terial from its own CPD or from the contiguous disc feature.Our models are in agreement with the previous mass rangesestimated for the two companions; further observations andmodelling of the disc and their atmospheres are needed tobetter constrain their masses.The estimated accretion rates of the companions are ofthe order of ∼ − M (cid:12) · yr − (Haffert et al. 2019), thusradiation from accretion shocks near both planets are neg-ligible. From our results, accretion flux would only have anincidence in the modelled IR planet fluxes at distances ∼ , since accretion rates are expected to be higher due tothe scaling. This can be appreciated in the contrast curveof the three planet models in H -band: planets’ contrasts de-crease at these distances. The effect of the accretion shock’sradiation becomes negligible at (cid:38)
10 AU . HL Tau is one of the most extensively studied protoplane-tary discs, with several rings and gaps detected in the dustcontinuum (ALMA Partnership et al. 2015). It is a youngstellar object of ≤ at around
140 pc to us (Kenyonet al. 2008), with an estimated stellar mass of ∼ . M (cid:12) (Kenyon & Hartmann 1995; Close et al. 1997). Observationswere carried out using the LBTI L/M IR Camera (LMIR-cam, Skrutskie et al. 2010; Leisenring et al. 2012), usingonly one of the two primary mirrors of the LBT telescope.No point-sources were detected. For the normalisation of thesurface density, we took the inferred gas surface density fromCARMA observations (Kwon et al. 2011, 2015) at a fiducial MNRAS000
140 pc to us (Kenyonet al. 2008), with an estimated stellar mass of ∼ . M (cid:12) (Kenyon & Hartmann 1995; Close et al. 1997). Observationswere carried out using the LBTI L/M IR Camera (LMIR-cam, Skrutskie et al. 2010; Leisenring et al. 2012), usingonly one of the two primary mirrors of the LBT telescope.No point-sources were detected. For the normalisation of thesurface density, we took the inferred gas surface density fromCARMA observations (Kwon et al. 2011, 2015) at a fiducial MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations Figure 11.
Application of our model to planets with . , . and . M J embedded in the PDS disc. The contrast curvesshown for H -, K - and L -bands were obtained considering stellarmagnitudes of H = . , K = . and L = . (Cutriet al. 2003; Cutri & et al. 2014). The two planetary companions(Keppler et al. 2018; Mesa et al. 2019; Haffert et al. 2019) areshown as black and grey crosses, with the corresponding uncer-tainties. distance of
40 AU , Σ =
34 g · cm − (a factor × . comparedto the simulated disc).In Figure 12 we show the contrast limit of the LBTIobservation in L -band as a function of the angular separa-tion to the central star, together with the derived contrastof the re-scaled models for planets with . , . and . M J . In Appendix B, contrast curves in J -, H - and K -bandsare included for completeness. In L -band, a high extinction Figure 12.
Contract curves in L -band for planets embeddedin HL Tau, including the σ detection limit of the observationfrom Testi et al. (2015). The observations were performed us-ing LBTI/LMIRcam. The contrast curves are for planet massesof . , . and . M J . The considered apparent magnitude ofthe central star was L = .
23 mag (Testi et al. 2015). The colouredregions accounts for the uncertainty in the planet contrast. Thegrey vertical area is delimited by the D and D rings detectedin dust continuum (ALMA Partnership et al. 2015). is predicted for . and . M J along the entire disc, es-pecially at distances (cid:46)
60 AU ; at that distance, A L valuesare .
27 mag ( A V = .
29 mag ) and .
25 mag ( A V = . ) for these planets respectively. For planets outer in thedisc, the extinction contribution is smaller but still signifi-cant: .
02 mag ( A V = .
29 mag ) for . M J , and .
88 mag ( A V = .
13 mag ) for . M J at
120 AU . These planets arenot massive enough to clear the gap efficiently. On the otherhand, for the most massive planet ( . M J ), extinction isnegligible at any distance.Six gaps were observed in the ALMA continuum ob-servation; following the example as in Testi et al. (2015),within the gap delimited by D5 and D6 rings (marked asgrey vertical line) the contrast limit of the instrument doesnot allow us to constrain the mass of the companion thatcould be responsible for the gap. Nevertheless, from the in-ferred contrast curves, extinction would have an incidence ina hypothetical point-source detection only for planet masses (cid:46) . M J . Observations using the Keck/NIRC2 vortex coronagraphwere performed by Ruane et al. (2017) searching for point-sources in the TW Hya disc. This system is the closest knownprotoplanetary disc to us ( . , Gaia Collaboration et al.2018), with a central star of . - . M (cid:12) (Andrews et al. 2012;Herczeg & Hillenbrand 2014) relatively old ( -
10 Myr , Ruaneet al. 2017), and with an estimated total disc mass of . M (cid:12) (Bergin et al. 2013). The surface gas density has been mod-elled from ALMA line emission observations (Kama et al.2016; Trapman et al. 2017). From their unperturbed mod-els, we used a fiducial surface density of . · cm − at for the re-scaling (a surface density factor × . of thesimulated disc). The instrument allows for IR high-contrastimaging in L -band, using angular differential imaging (ADI) MNRAS , 1–18 (2020) Sanchis et al.
Figure 13.
Contract curves in L -band for planets embed-ded in TW Hya, including significance detection limits ofKeck/NIR observations (Ruane et al. 2017). The contrast limitsare shown for angular differential imaging (ADI) and referencestar differential imaging (RDI). The contrast curves shown arefor planet masses of . , . and . M J . The apparent mag-nitude of the central star is L = .
01 mag , taken from the W1band in the WISE catalogue (Wright et al. 2010). The colouredregions for each planet model are delimited by the estimated ages( -
10 Myr , Ruane et al. 2017). The grey vertical lines account forthe gaps observed in Andrews et al. (2016) and van Boekel et al.(2017). and reference star differential imaging (RDI). In Figure 13,the detection limits for ADI and RDI are shown togetherwith the expected contrast of planets with . , . and . M J . In this observation, RDI allows for detections ofpoint-sources at distances as low as ∼ . At this distancethe accretion flux overcomes the intrinsic flux for a planetof (cid:38) . M J , thus the contrast decreases compared to thenon-accreting case. The contrast curves of these planets in J -, H - and K -bands are shown in Appendix B.Different observations of TW Hya confirmed severalgaps in the disc. Andrews et al. (2016) detected three darkannuli at , and
47 AU distances to the host star (dis-tances corrected with newest Gaia parallax). An unresolvedgap in the inner disc was also seen from µ m contin-uum emission using ALMA. Scattered light using SPHEREdetected three gaps in the polarised intensity distribution,at (cid:46) , , and
88 AU . In the figure we show the gaps inthe outer disc ( ∼ , , and
87 AU ). No point-sourceswere detected in the Keck/NIRC2 observations. Ruane et al.(2017) set upper limits for planets located at these gaps ( . - . M J , . - . M J , . - . M J , and . - . M J from innerto outer distances). Analogously, we can infer upper limitsof the planets interpolating our results, since the contrastcurves lay between our models. Using the models for planets, the upper limits for these gaps would be . M J , . M J , . M J and . M J . Considering an age of
10 Myr ,the upper limits are marginally higher: . M J , . M J , . M J and . M J . Taking into account the extinction due tothe disc above the planet increases the estimated upper lim-its compared to the previous work, which did not considerthis effect. This indicates the importance of extinction whenlooking for protoplanets with direct imaging methods. Figure 14.
Contract curves in L -band for planets embedded inHD , including the σ detection limits of the observationfrom Guidi et al. (2018). The observations were performed usingKeck/NIRC2. The contrast of planets with . , . and . M J are shown. The apparent magnitude of the central star is L = . , inferred from the W1 band in the WISE catalogue (Wrightet al. 2010). The grey vertical lines account for the gaps observedin Isella et al. (2016, 2018) . In Guidi et al. (2018), the HD disc was studied inthe L -band using the same instrument (Keck/NIRC2 vor-tex coronagraph). The scattered polarised emission in the J -band was also studied with the Gemini Planet Imager inMonnier et al. (2017), detecting a ring with an offset thatcan be explained by an inclined flared disc. This system hasa central star of . M (cid:12) (Natta et al. 2004) and estimated ageof ∼ (Montesinos et al. 2009). Observations in the dustcontinuum using ALMA (Isella et al. 2016) confirmed the ex-istence of three gaps at distances of ∼ , ∼ and ∼
136 AU (corrected with the new Gaia distance of . ). Kine-matical analysis of gas observations suggested the presenceof two planets at the second and third gaps (Teague et al.2018). In Pinte et al. (2018), HD models showed that a thirdplanet is expected further out. The estimated masses of thethree potential planets are M J (at
83 AU ), . M J (at ) and ≈ M J at ( ≈
260 AU ). The new DSHARP/ALMAobservations confirmed an additional gap at ∼
10 AU (Isellaet al. 2018); assuming that this gap is caused by a planet,Zhang et al. (2018) estimated a planet mass between . and . M J from 2D HD simulations.The L -band high-contrast imaging (Guidi et al. 2018)detected a point-like source at a distance of . with . σ significance. None of the observations in L - or J -bandfound any point-sources at the gaps observed in the con-tinuum. Our models allow to set upper limits for planetsat the location of the gaps. We used a fiducial surface den-sity of Σ = . · cm − at
40 AU (from Isella et al. 2016),corresponding to a factor × . of the simulated disc, to ob-tain contrast curves for . , . and . M J ( L -band inFigure 14, and J -, H - and K -bands in Figures B4 in the Ap-pendix B). The innermost gap in Figure 14 is within themasked region in the Keck/NIR2 observations, thus a massupper limit can not be inferred. At the second gap, the modelfor our most massive planet lays slightly below the detec- MNRAS000
40 AU (from Isella et al. 2016),corresponding to a factor × . of the simulated disc, to ob-tain contrast curves for . , . and . M J ( L -band inFigure 14, and J -, H - and K -bands in Figures B4 in the Ap-pendix B). The innermost gap in Figure 14 is within themasked region in the Keck/NIR2 observations, thus a massupper limit can not be inferred. At the second gap, the modelfor our most massive planet lays slightly below the detec- MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations tion limit of the observation. A rough extrapolation wouldyield an upper-limit of . M J , slightly below the range pro-vided by Guidi et al. (2018) ( - M J in that work). Forthe third and fourth gaps, we obtain upper limits of . M J and . M J from interpolating our models. These values areslightly higher than the upper limits inferred in Guidi et al.(2018) ( . - . M J , and . - M J respectively). Taking intoaccount extinction on the contrast of the planets increasethe inferred upper limits of the non-detected planets. In ev-ery gap, extinction does have an important effect for planetswith masses lower than the inferred upper limits. Comparedto the estimates of Teague et al. (2018) and Pinte et al.(2018) from indirect analysis, our inferred upper limits aresignificantly higher; consequently a direct detection of thesecompanions would only be possible improving the detectionlimit to much higher contrast. In this work we studied the effect of extinction for directimaging of young planets embedded in protoplanetary discs.A set of HD simulations were performed to reproduce planet-disc interaction at high resolution for several planet masses.Column densities and extinction coefficients were derivedin order to model the planet predicted magnitudes in J -, H -, K -, L -, M -, N -bands. Exploiting properties of locallyisothermal discs, we applied the models to planets embed-ded in CQ Tau, PDS and HL Tau protoplanetary discs,and inferred upper-limits for planets at the gaps observedin TW Hya and HD . The most important results ofthis work are: • For the simulated planets at . , the M J clears itssurrounding material very effectively. The resulting columndensity is extremely low, and, as consequence, extinction isnot significant in any band. The and M J planets arecompletely hidden by the disc at (cid:46) µ m wavelengths (withrespective extinction coefficients of > , >
15 mag ), whileat wavelengths between to µ m their corresponding co-efficients are reduced, below and . In the N -band,extinction is higher compared to L - and N -bands due to thesilicate feature in the assumed ISM dust opacities. • Jupiter-like planets embedded in discs with very lowunperturbed surface densities (of the order of (cid:46) · cm − )have very low extinction coefficients in IR at any distanceconsidered. In CQ Tau, only planets with (cid:46) M J are affectedby extinction in J - and H -bands at distances (cid:46)
20 AU . InPDS , extinction has an incidence only for the least mas-sive planet model at distances <
50 AU , more significant atshorter wavelengths. • In more dense discs like HL Tau and HD , di-rect detection of companions is unlikely in J -, H -, K -, and N -bands due to the extinction effects. Only the most mas-sive planet from our models would be detectable, since theextinction is negligible. • We inferred upper limits of the gaps in TW Hya andHD , slightly higher than previous work due to theeffect of extinction. This points out the importance of ex-tinction from the disc material in high-contrast imaging ofprotoplanetary discs. • Radiation from accretion shocks onto the planet hasbeen considered in our models. It can have an important effect on the total planet emission for accretion rates of theorder of ∼ − M (cid:12) / yr ; these high rates occur at distances (cid:46)
10 AU in our models.The scarcity of detections so far might suggest differentscenarios: giant planet formation further out in the disc israre, or perhaps planets formed at these early stages are stillnot massive enough ( (cid:46) M J ) to be detected with currentinstrumentation. ACKNOWLEDGEMENTS
REFERENCES
ALMA Partnership et al., 2015, ApJ, 808, L3Andrews S. M., 2015, PASP, 127, 961Andrews S. M., et al., 2012, ApJ, 744, 162Andrews S. M., et al., 2016, ApJ, 820, L40Andrews S. M., et al., 2018, ApJ, 869, L41Ansdell M., Williams J. P., Manara C. F., Miotello A., FacchiniS., van der Marel N., Testi L., van Dishoeck E. F., 2017, AJ,153, 240Beccari G., et al., 2010, ApJ, 720, 1108Bergin E. A., et al., 2013, Nature, 493, 644Brown J. M., Blake G. A., Qi C., Dullemond C. P., Wilner D. J.,2008, ApJ, 675, L109Bryden G., Chen X., Lin D. N. C., Nelson R. P., PapaloizouJ. C. B., 1999, ApJ, 514, 344Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245Chapillon E., Guilloteau S., Dutrey A., Pi´etu V., 2008, A&A, 488,565Chiar J. E., Tielens A. G. G. M., 2006, ApJ, 637, 774Close L. M., Roddier F., J. Northcott M., Roddier C., Elon GravesJ., 1997, ApJ, 478, 766Cutri R. M., et al. 2014, VizieR Online Data Catalog, 2328Cutri R. M., et al., 2003, VizieR Online Data Catalog, 2246D’Angelo G., Kley W., Henning T., 2003, ApJ, 586, 540Dipierro G., Laibe G., 2017, MNRAS, 469, 1932Dong R., Fung J., 2017, ApJ, 835, 146Dong R., Zhu Z., Whitney B., 2015a, ApJ, 809, 93Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015b, ApJ, 809,L5Duffell P. C., 2015, ApJ, 806, 182D¨urmann C., Kley W., 2015, A&A, 574, A52MNRAS , 1–18 (2020) Sanchis et al.
Ercolano B., Pascucci I., 2017, Royal Society Open Science, 4,170114Frank J., King A. R., Raine D. J., 1985, Accretion power in as-trophysicsFung J., Dong R., 2015, ApJ, 815, L21Gaia Collaboration et al., 2018, A&A, 616, A1Guidi G., et al., 2018, MNRAS, 479, 1505G¨uver T., ¨Ozel F., 2009, MNRAS, 400, 2050Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinch-mann J., Girard J. H., Keller C. U., Bacon R., 2019, NatureAstronomy, 3, 749Hartmann L., 1998, Accretion Processes in Star FormationHartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135Hayashi C., 1981, Progress of Theoretical Physics Supplement,70, 35Herczeg G. J., Hillenbrand L. A., 2014, ApJ, 786, 97Hern´andez J., Hartmann L., Calvet N., Jeffries R. D., GutermuthR., Muzerolle J., Stauffer J., 2008, ApJ, 686, 1195Isella A., et al., 2016, Physical Review Letters, 117, 251101Isella A., et al., 2018, ApJ, 869, L49Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S.,P´erez L., 2019, ApJ, 879, L25Jang-Condell H., 2009, ApJ, 700, 820Jang-Condell H., 2017, ApJ, 835, 12Jang-Condell H., Turner N. J., 2013, ApJ, 772, 34Juh´asz A., Rosotti G. P., 2018, MNRAS, 474, L32Juh´asz A., Benisty M., Pohl A., Dullemond C. P., Dominik C.,Paardekooper S. J., 2015, MNRAS, 451, 1147Kama M., et al., 2016, A&A, 592, A83Kenyon S. J., Hartmann L., 1995, ApJS, 101, 117Kenyon S. J., G´omez M., Whitney B. A., 2008, Low Mass StarFormation in the Taurus-Auriga Clouds. p. 405Keppler M., et al., 2018, preprint, ( arXiv:1806.11568 )Klahr H., Kley W., 2006, A&A, 445, 747Kley W., 1999, MNRAS, 303, 696Kwon W., Looney L. W., Mundy L. G., 2011, ApJ, 741, 3Kwon W., Looney L. W., Mundy L. G., Welch W. J., 2015, ApJ,808, 102Leisenring J. M., et al., 2012, in Ground-based and Air-borne Instrumentation for Astronomy IV. p. 84464F,doi:10.1117/12.924814Manara C. F., Morbidelli A., Guillot T., 2018, A&A, 618, L3Mendigut´ıa I., Calvet N., Montesinos B., Mora A., Muzerolle J.,Eiroa C., Oudmaijer R. D., Mer´ın B., 2011, A&A, 535, A99Meru F., Rosotti G. P., Booth R. A., Nazari P., Clarke C. J.,2019, MNRAS, 482, 3678Mesa D., et al., 2019, A&A, 632, A25Mignone A., Tzeferacos P., Zanni C., Tesileanu O., MatsakosT., Bodo G., 2010, PLUTO: A Code for Flows in Multi-ple Spatial Dimensions, Astrophysics Source Code Library(ascl:1010.045)Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P.,Bodo G., 2012, ApJS, 198, 7Miotello A., et al., 2017, A&A, 599, A113Monnier J. D., et al., 2017, ApJ, 838, 20Montesinos B., Eiroa C., Mora A., Mer´ın B., 2009, A&A, 495, 901M¨uller A., et al., 2018, A&A, 617, L2Natta A., Testi L., Neri R., Shepherd D. S., Wilner D. J., 2004,A&A, 416, 179Nazari P., Booth R. A., Clarke C. J., Rosotti G. P., Tazzari M.,Juhasz A., Meru F., 2019, MNRAS, 485, 5914Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435,2610Ormel C. W., Paszun D., Dominik C., Tielens A. G. G. M., 2009,A&A, 502, 845Ormel C. W., Min M., Tielens A. G. G. M., Dominik C., PaszunD., 2011, A&A, 532, A43Ormel C. W., Kuiper R., Shi J.-M., 2015a, MNRAS, 446, 1026 Ormel C. W., Shi J.-M., Kuiper R., 2015b, MNRAS, 447, 3512Paardekooper S.-J., 2017, MNRAS, 469, 4306Paardekooper S. J., Mellema G., 2006, A&A, 453, 1129Pi´etu V., Dutrey A., Guilloteau S., Chapillon E., Pety J., 2006,A&A, 460, L43Pinte C., et al., 2018, ApJ, 860, L13Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., RoushT., Fong W., 1994, ApJ, 421, 615Quanz S. P., Amara A., Meyer M. R., Girard J. H., KenworthyM. A., Kasper M., 2015, ApJ, 807, 64Rameau J., Chauvin G., Lagrange A.-M., Maire A.-L., BoccalettiA., Bonnefoy M., 2015, A&A, 581, A80Reggiani M., et al., 2014, ApJ, 792, L23Reggiani M., et al., 2018, A&A, 611, A74Rice W. K. M., Armitage P. J., Wood K., Lodato G., 2006, MN-RAS, 373, 1619Rosotti G. P., Juhasz A., Booth R. A., Clarke C. J., 2016, MN-RAS, 459, 2790Ruane G., et al., 2017, AJ, 154, 73Scicluna P., Rosotti G., Dale J. E., Testi L., 2014, A&A, 566, L3Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Skrutskie M. F., et al., 2010, in Ground-based and Air-borne Instrumentation for Astronomy III. p. 77353H,doi:10.1117/12.857724Spiegel D. S., Burrows A., 2012, ApJ, 745, 174Szul´agyi J., Morbidelli A., Crida A., Masset F., 2014, ApJ, 782,65Szul´agyi J., Masset F., Lega E., Crida A., Morbidelli A., GuillotT., 2016, MNRAS, 460, 2853Tanigawa T., Watanabe S.-i., 2002, ApJ, 580, 506Tazzari M., et al., 2016, A&A, 588, A53Tazzari M., et al., 2017, A&A, 606, A88Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-MackeyD., 2018, ApJ, 860, L12Testi L., Natta A., Shepherd D. S., Wilner D. J., 2003, A&A, 403,323Testi L., et al., 2015, ApJ, 812, L38Trapman L., Miotello A., Kama M., van Dishoeck E. F., BrudererS., 2017, A&A, 605, A69Trotta F., Testi L., Natta A., Isella A., Ricci L., 2013, A&A, 558,A64Ubeira Gabellini M. G., et al., 2019, MNRAS, 486, 4638Wolf S., D’Angelo G., 2005, ApJ, 619, 1114Wright E. L., et al., 2010, AJ, 140, 1868Zhang S., et al., 2018, ApJ, 869, L47Zhu Z., 2015, ApJ, 799, 16van Boekel R., et al., 2017, ApJ, 837, 132
APPENDIX A: MAGNITUDES FOR PLANETSAT 10, 20, 50 AND 100 AU
The predicted magnitudes of the modelled planets at variousdistances to the central star are included in Tables A1, A2,A3 and A4 for completeness.
APPENDIX B: CONTRAST OF PLANETSEMBEDDED IN CQ TAU, HL TAU, TW HYAAND HD 163296
The remaining contrast curves of planets as a function ofthe distance to the host star in J -, H - and K -bands of all thestudied systems: . , . and . M J planets in CQ Taudisc (Figure B1), . , . and . M J planets in HL Tau MNRAS000
The remaining contrast curves of planets as a function ofthe distance to the host star in J -, H - and K -bands of all thestudied systems: . , . and . M J planets in CQ Taudisc (Figure B1), . , . and . M J planets in HL Tau MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations Table A1.
Expected magnitudes for planets at
10 AU to the host star, with masses: M J –viscous and inviscid scenarios– M J and M J . Mag pl is the total magnitude of the planet, including accretion flux; A band is the extinction coefficient in each band, Mag expected is thepredicted magnitude of the planet considering extinction due to disc material.
Hot-start planet
Cold-start planetJ H K L M N J H K L M N
Mag pl .
13 14 .
95 13 .
31 12 .
65 11 .
40 10 .
02 15 .
92 16 .
02 15 .
95 14 .
94 12 .
97 12 . M J A band .
58 41 .
30 25 .
99 14 .
11 11 .
05 19 .
91 64 .
58 41 .
30 25 .
99 14 .
11 11 .
05 19 . expected .
70 56 .
25 39 .
30 26 .
76 22 .
45 29 .
93 80 .
49 57 .
32 41 .
94 29 .
05 24 .
03 31 . pl .
55 12 .
76 11 .
66 11 .
01 10 .
59 9 .
34 14 .
93 15 .
01 15 .
00 14 .
12 12 .
64 11 . M J A band .
62 11 .
91 7 .
49 4 .
13 3 .
23 5 .
82 18 .
62 11 .
91 7 .
49 4 .
13 3 .
23 5 . expected .
17 24 .
67 19 .
16 15 .
14 13 .
82 15 .
16 33 .
55 26 .
92 22 .
50 18 .
25 15 .
88 17 . pl .
25 10 .
28 9 .
52 9 .
04 9 .
33 8 .
35 14 .
80 14 .
84 14 .
89 13 .
52 12 .
39 11 . M J A band .
26 0 .
16 0 .
10 0 .
06 0 .
04 0 .
08 0 .
26 0 .
16 0 .
10 0 .
06 0 .
04 0 . expected .
51 10 .
44 9 .
62 9 .
10 9 .
38 8 .
43 15 .
05 15 .
00 14 .
99 13 .
58 12 .
44 11 . pl .
53 15 .
25 13 .
36 12 .
68 11 .
41 10 .
03 17 .
51 17 .
69 17 .
27 15 .
33 13 .
02 12 . M inviscid A band .
86 17 .
82 11 .
21 6 .
18 4 .
84 8 .
71 27 .
86 17 .
82 11 .
21 6 .
18 4 .
84 8 . expected .
39 33 .
07 24 .
58 18 .
86 16 .
25 18 .
74 45 .
37 35 .
51 28 .
49 21 .
51 17 .
86 20 . Table A2.
Absolute magnitudes for planets at
20 AU to the host star, with masses: M J –viscous and inviscid scenarios– M J and M J . Mag pl is the total magnitude of the planet, including accretion flux; A band is the extinction coefficient in each band, Mag expected is thepredicted magnitude of the planet considering extinction due to disc material.
Hot-start planet
Cold-start planetJ H K L M N J H K L M N
Mag pl .
55 15 .
27 13 .
37 12 .
68 11 .
41 10 .
03 17 .
77 17 .
98 17 .
45 15 .
36 13 .
02 12 . M J A band .
66 29 .
20 18 .
38 9 .
98 7 .
82 14 .
08 45 .
66 29 .
20 18 .
38 9 .
98 7 .
82 14 . expected .
22 44 .
47 31 .
74 22 .
66 19 .
22 24 .
10 63 .
43 47 .
18 35 .
83 25 .
34 20 .
84 26 . pl .
76 12 .
85 11 .
69 11 .
03 10 .
60 9 .
34 16 .
86 16 .
91 16 .
70 14 .
64 12 .
74 11 . M J A band .
17 8 .
42 5 .
30 2 .
92 2 .
29 4 .
12 13 .
17 8 .
42 5 .
30 2 .
92 2 .
29 4 . expected .
93 21 .
28 16 .
99 13 .
95 12 .
89 13 .
46 30 .
03 25 .
33 22 .
00 17 .
56 15 .
02 15 . pl .
28 10 .
29 9 .
52 9 .
05 9 .
34 8 .
35 16 .
61 16 .
50 16 .
55 14 .
16 12 .
55 11 . M J A band .
18 0 .
12 0 .
07 0 .
04 0 .
03 0 .
06 0 .
18 0 .
12 0 .
07 0 .
04 0 .
03 0 . expected .
46 10 .
40 9 .
59 9 .
09 9 .
37 8 .
41 16 .
79 16 .
61 16 .
62 14 .
20 12 .
58 11 . pl .
62 15 .
32 13 .
38 12 .
69 11 .
41 10 .
03 18 .
59 18 .
94 17 .
92 15 .
43 13 .
03 12 . M inviscid A band .
70 12 .
60 7 .
93 4 .
37 3 .
42 6 .
16 19 .
70 12 .
60 7 .
93 4 .
37 3 .
42 6 . expected .
32 27 .
92 21 .
30 17 .
06 14 .
83 16 .
19 38 .
29 31 .
54 25 .
84 19 .
79 16 .
45 18 . Table A3.
Expected magnitudes for planets at
50 AU to the host star, with masses: M J –viscous and inviscid scenarios– M J and M J . Mag pl is the total magnitude of the planet, including accretion flux; A band is the extinction coefficient in each band, Mag expected is thepredicted magnitude of the planet considering extinction due to disc material.
Hot-start planet
Cold-start planetJ H K L M N J H K L M N
Mag pl .
63 15 .
32 13 .
38 12 .
69 11 .
41 10 .
03 18 .
77 19 .
18 18 .
00 15 .
43 13 .
03 12 . M J A band .
88 18 .
47 11 .
62 6 .
31 4 .
94 8 .
90 28 .
88 18 .
47 11 .
62 6 .
31 4 .
94 8 . expected .
51 33 .
79 25 .
00 19 .
00 16 .
35 18 .
93 47 .
65 37 .
65 29 .
62 21 .
75 17 .
97 20 . pl .
80 12 .
87 11 .
70 11 .
03 10 .
60 9 .
34 18 .
03 18 .
00 17 .
46 14 .
73 12 .
75 11 . M J A band .
33 5 .
33 3 .
35 1 .
85 1 .
45 2 .
60 8 .
33 5 .
33 3 .
35 1 .
85 1 .
45 2 . expected .
12 18 .
19 15 .
05 12 .
88 12 .
05 11 .
94 26 .
35 23 .
33 20 .
81 16 .
58 14 .
20 14 . pl .
28 10 .
29 9 .
52 9 .
05 9 .
34 8 .
35 17 .
53 17 .
21 17 .
27 14 .
28 12 .
57 11 . M J A band .
11 0 .
07 0 .
05 0 .
03 0 .
02 0 .
04 0 .
11 0 .
07 0 .
05 0 .
03 0 .
02 0 . expected .
39 10 .
36 9 .
57 9 .
08 9 .
36 8 .
39 17 .
65 17 .
28 17 .
31 14 .
31 12 .
59 11 . pl .
64 15 .
33 13 .
38 12 .
69 11 .
41 10 .
03 18 .
87 19 .
32 18 .
04 15 .
44 13 .
03 12 . M inviscid A band .
46 7 .
97 5 .
01 2 .
76 2 .
16 3 .
90 12 .
46 7 .
97 5 .
01 2 .
76 2 .
16 3 . expected .
10 23 .
30 18 .
39 15 .
45 13 .
57 13 .
93 31 .
32 27 .
29 23 .
05 18 .
20 15 .
19 15 . MNRAS , 1–18 (2020) Sanchis et al.
Table A4.
Expected magnitudes for planets at
100 AU to the host star, with masses: M J –viscous and inviscid scenarios– M J and M J . Mag pl is the total magnitude of the planet, including accretion flux; A band is the extinction coefficient in each band, Mag expected is thepredicted magnitude of the planet considering extinction due to disc material.
Hot-start planet
Cold-start planetJ H K L M N J H K L M N
Mag pl .
63 15 .
33 13 .
38 12 .
69 11 .
41 10 .
03 18 .
87 19 .
33 18 .
04 15 .
44 13 .
03 12 . M J A band .
42 13 .
06 8 .
22 4 .
46 3 .
50 6 .
30 20 .
42 13 .
06 8 .
22 4 .
46 3 .
50 6 . expected .
06 28 .
39 21 .
60 17 .
15 14 .
90 16 .
32 39 .
29 32 .
39 26 .
26 19 .
90 16 .
52 18 . pl .
80 12 .
87 11 .
70 11 .
03 10 .
60 9 .
34 18 .
16 18 .
12 17 .
53 14 .
74 12 .
75 11 . M J A band .
89 3 .
77 2 .
37 1 .
31 1 .
02 1 .
84 5 .
89 3 .
77 2 .
37 1 .
31 1 .
02 1 . expected .
69 16 .
64 14 .
07 12 .
34 11 .
62 11 .
18 24 .
05 21 .
89 19 .
90 16 .
04 13 .
77 13 . pl .
28 10 .
29 9 .
52 9 .
05 9 .
34 8 .
35 17 .
63 17 .
27 17 .
33 14 .
29 12 .
57 11 . M J A band .
08 0 .
05 0 .
03 0 .
02 0 .
01 0 .
03 0 .
08 0 .
05 0 .
03 0 .
02 0 .
01 0 . expected .
36 10 .
34 9 .
55 9 .
07 9 .
35 8 .
38 17 .
71 17 .
32 17 .
36 14 .
31 12 .
58 11 . pl .
64 15 .
33 13 .
38 12 .
69 11 .
41 10 .
03 18 .
89 19 .
34 18 .
05 15 .
44 13 .
03 12 . M inviscid A band .
81 5 .
63 3 .
55 1 .
95 1 .
53 2 .
76 8 .
81 5 .
63 3 .
55 1 .
95 1 .
53 2 . expected .
45 20 .
96 16 .
92 14 .
64 12 .
94 12 .
79 27 .
70 24 .
98 21 .
59 17 .
39 14 .
56 14 . (Figure B2), . , . and . M J planets in TW Hya (Fig-ure B3), and . , . and . M J planets in HD (Figure B4). This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000
56 14 . (Figure B2), . , . and . M J planets in TW Hya (Fig-ure B3), and . , . and . M J planets in HD (Figure B4). This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000 , 1–18 (2020) etectability of protoplanets from HD simulations Figure B1.
Contrast in J -, H - and K -bands of planets embeddedin the CQ Tau disc as a function of distance to the central star .Results are fot planets with . , . and . M J . Figure B2.
Application of our model to planets with . , . and . M J masses embedded in the HL Tau disc. Contrast isshown for different distances to the host star in J -, H - and K -bands.MNRAS , 1–18 (2020) Sanchis et al.
Figure B3.
Application of our model to planets with . , . and . M J masses embedded in the TW Hya disc. Contrastis shown for different distances to the host star in J -, H - and K -bands. Figure B4.
Application of our model to planets with . , . and . M J masses embedded in the HD disc. Contrastis shown for different distances to the host star in J -, H - and K -bands. MNRAS000