Keplerian motion of a compact source orbiting the inner disc of PDS 70: a third protoplanet in resonance with b and c?

ABSTRACT

The disc around PDS 70 hosts two directly imaged protoplanets in a gap. Previous VLT/SPHERE and recent JWST/NIRCam observations have hinted at the presence of a third compact source in the same gap at |$\sim$|13 au, interior to the orbit of PDS 70 b. We reduce seven published and one unpublished VLT/SPHERE data sets in YJH and K bands, as well as an archival VLT/NaCo data set in L’ band, and an archival VLT/SINFONI data set in H + K band. We combine angular-, spectral-, and reference star differential imaging to search for protoplanet candidates. We recover the compact source in all epochs, consistent with the JWST detection, moving on an arc that can be fit by Keplerian motion of a protoplanet that could be in resonance with PDS 70 b and c. We find that the spectral slope is overall consistent with the unresolved star and inner disc emission at 0.95–1.65 |$\mu\mathrm{m}$|⁠, which suggests a dust-scattering-dominated spectrum. An excess beyond 2.3 |$\mu\mathrm{m}$| could be thermal emission from either a protoplanet or heated circumplanetary dust, variability, or inner disc contamination, and requires confirmation. While we currently cannot rule out a moving inner disc feature or a dust clump associated with an unseen planet, the data support the hypothesis of a third protoplanet in this remarkable system.

planet–disc interactions, protoplanetary discs, stars: individual: PDS 70, infrared: planetary systems

1 INTRODUCTION

Mean-motion resonance (MMR) capture is expected to be a key mechanism that sets the orbital architecture of protoplanets during the disc phase and halts migration towards the stellar potential well (Krijt et al. 2023; Weiss et al. 2023). MMR is observed in the Galilean satellites Io, Europa, and Ganymede around Jupiter (Laplace 1799) and would explain the coplanar 4:2:1 resonance of the Gliese-876 exoplanets (Nelson et al. 2016; Cimerman, Kley & Kuiper 2018), of which two of the planets are several times the mass of Jupiter. Imaging of the young (⁠|$\sim$|30–60 Myr) system HR 8799 revealed four giant planets consistent with a coplanar 8:4:2:1 Laplace resonance for long-term stability (Marois et al. 2008, 2010; Fabrycky & Murray-Clay 2010; Goździewski & Migaszewski 2020). However, our understanding of how and when protoplanets enter such a configuration is limited due to the small number of confirmed multiprotoplanet systems.

PDS 70 (V1032 Centauri) is a |$\sim$|5Myr, 0.87–1|$M_\odot$| T Tauri star (Hashimoto et al. 2012; Long et al. 2018; Müller et al. 2018; Wang et al. 2021) at a distance of 112.39|$\pm$|0.24 pc (Gaia Collaboration 2023) surrounded by a protoplanetary disc. The system hosts two directly imaged protoplanets (hereafter PDS 70 b and c) with infrared emission, H|$\alpha$| emission, and a sub-mm continuum counterpart (Keppler et al. 2018; Müller et al. 2018; Haffert et al. 2019; Benisty et al. 2021). PDS 70 b and c have a deprojected separation of |$\sim$|20 and 34 au from the star, with orbital fits of b preferring non-zero eccentricity (Wang et al. 2021). Dedicated hydrodynamical simulations found that the protoplanets lock into a 2:1 MMR (Bae et al. 2019), and migrate to the same configuration even if initially placed on wider orbits (Toci et al. 2020).

Mesa et al. (2019) first identified a candidate third protoplanet interior to the orbit of b at |$\sim$|13 au, with the Spectro-Polarimetic High contrast imager for Exo-planets REsearch (SPHERE; Beuzit et al. 2008) instrument at the VLT. Their ‘point-like feature’ (PLF) had a spectrum similar to the stellar spectrum in the range 0.95–1.65 |$\mu\mathrm{m}$| and it was marginally detected in H and K bands. Given the uncertain radial extent of the inner disc around PDS 70 A, they concluded that the PLF is likely forward scattering from disc material, although an embedded planet was not ruled out. Recent JWST/NIRCam observations with the F187N (1.87 |$\mu\mathrm{m}$|⁠) filter in Christiaens et al. (2024) reported a bright signal near the previously reported PLF and was robust to the subtraction of an optimal disc model. The measured astrometry is consistent with Keplerian motion of an object at |$\sim$|13.5 au in the plane of the disc. The feature was interpreted as either tracing a dusty circumplanetary envelope or disc, or a dust clump, and tentatively labelled PDS 70 d.

In this paper, we revisit archival data sets using novel and state-of-the-art high-contrast post-processing methods to constrain the nature and orbit of candidate protoplanet PDS 70 d. We present the high-contrast imaging data of PDS 70 and data reduction in Section 2. In Section 3, we present astrometry, orbital fitting, and the spectrum of candidate d. We discuss the potential planetary nature of the source in Section 4 and present our conclusions in Section 5.

2 OBSERVATIONS AND DATA REDUCTION

2.1 VLT/SPHERE IRDIFS observations

We processed six epochs of coronagraphic observations taken with the SPHERE instrument in IRDIFS mode, that is, simultaneous acquisition with the Infrared Differential Imaging Spectrometer (IRDIS) and the Integral Field Spectrograph (IFS) sub-instruments. These are summarized in Table 1. We reanalysed five observations that were presented in Mesa et al. (2019) and Wahhaj et al. (2024). The 2021 May 17 epoch is unpublished data from the SPHERE infrared survey for exoplanets (Desidera et al. 2021). Each sequence was taken in pupil stabilized mode, enabling the use of angular differential imaging (ADI; Marois et al. 2006a) with IRDIS and the IFS and spectral differential imaging (SDI; Sparks & Ford 2002) with the IFS. Three epochs from ESO program 107.22UJ.001 (PI: Z. Wahhaj, Wahhaj et al. 2024) observed a nearby PSF reference star, UCAC2 14412811, interspersed between every 10 min of integration time on PDS 70, enabling ‘star-hopping’ (Wahhaj et al. 2021) reference-star differential imaging (RDI).

Table 1.

Summary of high-contrast imaging observations of PDS 70 reduced and analysed in this work. FITS files are in the Data Availability statement.

Date	Strategy	Program ID	Instrument	Filter/	Plate scale|$^{\textit{a}}$|	Coronagraph	DIT	T|$_{\rm int}^{^{\textit{b}}}$|	|$\langle \epsilon \rangle ^{\textit{c}}$|	|$\Delta$|PA|$^{\textit{d}}$|
				Prism	(mas px|$^{-1}$|⁠)		(s)	(min)	(⁠|$^{\prime \prime }$|⁠)	(⁠|$^\circ$|⁠)
2012 Mar 31|$^{\textit{e}}$|	ADI	GS-2012A-C-3	NICI	L’	17.95|$\pm 0.01$|	–	0.76	33.57	–	99.37
2014 May 10	SADI	093.C-0526	SINFONI	K	12.5	–	1	116	1.13	99.8
2016 June 1	ADI	097.C-0206(A)	NaCo	L’	27.193|$\pm 0.09$|	–	0.2	62	0.5	85
2018 Feb 25	ADI	1100.C-0481(D)	IRDIS	K|$_1$|K|$_2$|	12.63|$\pm 0.009$|	N_ALC_YJH_S	96	112	0.4	91.86
2018 Feb 25	ASDI	1100.C-0481(D)	IFS	YJH	7.46|$\pm 0.02$|	N_ALC_YJH_S	96	136	0.4	91.86
2019 Mar 6	ASDI	1100.C-0481(L)	IFS	YJH	7.46|$\pm 0.02$|	N_ALC_Ks	96	76.8	0.39	56.4
2021 May 17	ASDI	1104.C-0416(E)	IFS	YJH	7.46|$\pm 0.02$|	N_ALC_YJH_S	96	76.8	0.46	56.38
2021 Aug 21	ASDI/RDI	107.22UJ.001	IFS	YJH	7.46|$\pm 0.02$|	N_ALC_YJH_S	32	27.2	0.75	14
2021 Sept 3	ASDI/RDI	107.22UJ.001	IFS	YJH	7.46|$\pm 0.02$|	N_ALC_YJH_S	32	21.3	0.54	5.1
2022 Feb 28	ADI/RDI	107.22UJ.001	IRDIS	Ks	12.265|$\pm 0.009$|	N_ALC_YJH_S	16	36.2	0.48	19.51
2022 Feb 28	ASDI/RDI	107.22UJ.001	IFS	YJH	7.46|$\pm 0.02$|	N_ALC_YJH_S	32	36.2	0.48	19.91
2023 Mar 8	Roll	1282	NIRCam	F187N	30.63|$\pm 0.06$|	–	0.35	8.3	–	5.0

Notes. |$^{\textit{a}}$|Maire et al. (2016) for IRDIS & IFS, Launhardt et al. (2020) for NaCo, and Keppler et al. (2018) for NICI. |$^{\textit{b}}$|Total integration time on PDS 70 after bad frame removal. |$^{\textit{c}}$|Average seeing at |$\lambda$|  = 500 nm over all PDS 70 integrations. |$^{\textit{d}}$|Parallactic angle variation after bad frame removal. |$^{\textit{e}}$|Appendix  A.

All SPHERE data were reduced with an in-house open-source pipeline, vcal-sphere¹ (Christiaens et al. 2023a), which uses esorex (v3.13.8) recipes for calibration and Vortex Image Processing (vip; Gomez Gonzalez et al. 2017; Christiaens et al. 2023b) functions for pre- and post-processing.

For IRDIS, the pipeline involved the same steps as presented in Hammond et al. (2023): (i) flat-fielding, (ii) principal component analysis (PCA)-based sky subtraction (Hunziker et al. 2018; Ren et al. 2023), (iii) bad pixel correction, (iv) centring of the coronagraphic images via satellite spots, (v) centring of the non-coronagraphic unsaturated PSF images using 2D Gaussian fits, (vi) anamorphism correction, (vii) bad-frame removal, (viii) final PSF creation, and (ix) post-processing with median-ADI (Marois et al. 2006b), PCA-ADI in full-frame (Amara & Quanz 2012; Soummer, Pueyo & Larkin 2012), annular PCA-ADI (Absil et al. 2013), and, when relevant, PCA-RDI. Each half of the detector was reduced independently for data sets taken in dual band imaging mode. As the 2018 February 25 ADI and 2022 February 28 RDI data sets had the longest integration time on target and excellent observing conditions, we consider only these two data sets to extract the photometry of candidate d in the |$K_1$||$K_2$| and |$K_s$| filters, respectively. The location of d near the inner working angle (IWA) of the coronagraph makes photometry extraction unreliable in data sets acquired in more challenging conditions.

Images obtained with PCA-ADI are affected by both self-subtraction and oversubtraction, which cause geometric artefacts to extended signals (see Juillard, Christiaens & Absil 2022 for the expected effect on the disc of PDS 70), while PCA-RDI reductions are still affected by oversubtraction. To alleviate these geometric artefacts, we also post-processed the 2018 February 25 and 2022 February 28 data sets using iterative PCA (IPCA; Pairet, Cantalloube & Jacques 2021)-ADI and IPCA-ARDI, respectively (Juillard, Christiaens & Absil 2023; Juillard et al. 2024). Given the stability of the 2018 data set and the close separation of d, the number of principal components (⁠|$n_{\rm pc}$|⁠) was set to 1 and the number of iterations (⁠|$n_{\rm it}$|⁠) to 10 for the IPCA-ADI reduction. For the 2022 data set, the science and reference stellar residual signals were not well correlated near the edge of the coronagraphic mask, hence we opted for an ARDI strategy with more principal components. Signal from d is detected with PCA-ARDI for |$n_{\rm pc}=2$|–16, with an optimal signal-to-noise ratio (SNR) achieved with 7 to 10 components. This reduction leads to significant self-subtraction and oversubtraction of the outer disc. For a faster convergence of IPCA-ARDI, we initialized its first iteration with the result of PCA-RDI obtained with data imputation (Ren 2023) instead of PCA-ARDI. As opposed to Ren (2023), we set the anchor mask to lie inside the cavity because in the case of PDS 70, disc signals extend up to the control ring. We selected an empty patch of signal in the cavity covering a wedge from PA = −30|$^{\circ }$| to PA = 15|$^{\circ }$| from 0.2 to 0.3 arcsec separation. We selected |$n_{\rm pc}=7$| for the IPCA-ARDI reduction, and the algorithm converged to a final image after |$n_{\rm it} = 126$| iterations. Finally, we also tested the mayonnaise² algorithm (Pairet et al. 2021) on the 2018 February 25 ADI data set, which attempts to separate point sources from extended circumstellar signals. The resulting image is shown in Fig. B1.

Calibration of IFS data included additional steps to the above. Master detector and instrument flats, spectra positions, and wavelength calibration were computed using esorex to produce calibrated spectral cubes of 39 channels. vip handled pre-processing in a similar fashion as for IRDIS. We leveraged the radial expansion of the satellite spots by |$\sim 14\lambda$|/D to determine the scaling factors used for SDI, and to identify channels with the highest SNR for the purposes of bad frame detection. We exploit both angular and SDI using PCA for improved subtraction of stellar speckles. We performed PCA in one step (i.e. PCA-ASDI) in order to reach the highest possible contrast in the absence of a dedicated reference star observation (e.g. Christiaens et al. 2021; Kiefer et al. 2021). A rotation threshold of 0.5|$\times$|FWHM at the separation of candidate d on each side of the considered frame was used to build the reference PSF library. We explore |$n_{\rm pc}=1\!-\!20$| and the frames in Fig. 1 are shown at |$n_{\rm pc}=10$|⁠; however we found that the speckles were sufficiently well subtracted after a few |$n_{\rm pc}$|⁠. Since images obtained with PCA-ASDI are also prone to geometrical artefacts (Christiaens et al. 2019a), we post-processed all data sets acquired in star-hopping mode with IPCA-RDI, applied channel per channel. We then median-combined the residual images, either band-wise (Y, J, and H) or over all spectral channels. For the data set acquired in the best conditions (2022 February 28), we also show the IPCA-RDI results obtained channel per channel in Fig. E1. For two out of the three data sets, the choice of |$n_{\rm pc}$| did not matter, as all tested values between 1 and 8 all led to very similar final images. For the 2021 August 21 data set, a minimum |$n_{\rm pc}$| value of 5 is required for a visual detection of candidate d, as the fluctuating observing conditions required a higher |$n_{\rm pc}$| value to capture most of the variance in the data. Most oversubtraction effects are then corrected within a few iterations. All IFS results obtained with IPCA-RDI are shown for |$n_{\rm pc} = 7$| and |$n_{\rm it} = 5$| in Fig. 1.

Figure 1.

9 yr of high-contrast imaging data of PDS 70. All images are cropped to a 1.2 arcsec |$\times$| 1.2 arcsec field of view centred on the star and in chronological order (left to right, top to bottom). The semitransparent mask over A is 160 mas wide, roughly the IWA of the SPHERE coronagraph (SINFONI, NaCo, and NIRCam are non-coronagraphic). The PSF-subtraction strategy is indicated in the top-left corner of the frame. The disc, b, and candidate d are detected in all frames, whereas c is clearest at |$\lambda >2~\mu\mathrm{m}$|⁠. North is up, east is to the left. Figure made using casa cube (https://github.com/cpinte/casa_cube).

Finally, we reanalysed two SPHERE data sets from 2015 May 3 and 2015 May 31. These data were taken in poor conditions and with the YJ prism instead of the wider YJH prism. The narrower spectral range reduces the effectiveness of SDI by excluding the H band. Since not even planet b was detected, we did not proceed with further analysis. As an additional test we processed both 2015 epochs with a second independent pipeline, charis-dep (Samland et al. 2022), with much the same result.

2.2 VLT/NaCo L’-band ADI observations

We made use of published non-coronagraphic VLT/NaCo L’-band (3.8 |$\mu\mathrm{m}$|⁠) observations from the ISPY survey (Launhardt et al. 2020) first presented in Keppler et al. (2018). These data were taken in pupil stabilized mode enabling the use of ADI. The data were acquired without a coronagraph, using a three-point dither pattern to subtract the thermal background. The sequence has a total integration time of 62 min on PDS 70 in good seeing (0.5 arcsec). We use PynPoint³ (Stolker et al. 2019) for the calibration and pre-processing of the data. Our pre-processing consists of (i) basic dark and flat calibrations, (ii) a simple thermal background subtraction using sky frames from other dither positions, (iii) a bad pixel interpolation, (iv) an alignment and centring of the star based on cross-correlations, and (v) a PCA-based bad frame rejection routine. After these pre-processing steps, we mean combined every five frames, resulting in the final cube of 1615 frames with 85|$^{\circ }$| of parallactic angle variation. We used annular PCA-ADI in a 2.7 arcsec subframe centred on the star and a rotation threshold of 0.5|$\times$|FWHM at the separation of PDS 70 d and explored |$n_{\rm pc}$|  = 1–50 with |$n_{\rm pc}$|  = 5 shown in Fig. 1.

2.3 VLT/SINFONI K-band SDI observations

We revisited VLT/SINFONI data from 2014 May 10 taken with the H  + K grating. The total integration time was 116 min (116 spectral cubes of 60s integration each), which resulted in 99.8° of parallactic angle variation in average observing conditions. Due to the 0.8 arcsec |$\times$| 0.8 arcsec field of view, a four-point dithering pattern was used where the star was placed in a different corner of the detector. Data reduction followed a similar approach as detailed in Christiaens et al. (2019a). None the less, compared to that work, we not only mean-combined quadruplets of spectral cubes obtained at different dither positions to enlarge the field of view of SINFONI, but also produced a different 4D (spectral + temporal) master cube considering a cropped field of view of 0.36 arcsec |$\times$| 0.36 arcsec. The latter was obtained by considering all individual spectral cubes where the star was at least 0.18 arcsec away from the edge of the field. Since the star drifted on the detector over the course of the observation, only 62 cubes complied to that condition. For both the dither-combined or the cropped master cubes, the 13 per cent least correlated cubes from the median spectral cube, corresponding to acquisition in the poorest conditions, were subsequently discarded. We then used the vip implementation of PCA-SADI to post-process the two different 4D data cubes, that is, PCA-SDI on each individual spectral cube of |$\sim$|2000 channels, followed by PCA-ADI on the cube of either 25 or 54 images resulting from PCA-SDI. The number of principal components was set to |$n_{\rm pc, S}=1$| or 2 for PCA-SDI, and |$n_{\rm pc, A}=1$|–10 was explored for PCA-ADI. Compared to Christiaens et al. (2019a), a smaller numerical mask was used to probe separations down to |$\sim$|62 mas (5 spaxels). While the full H + K cube of spectral images was used as PCA library for PCA-SDI, residual images after PCA model subtraction were only mean-combined at wavelengths longer than 2.3 |$\mu\mathrm{m}$| in the K band, where point-like signals benefit from the largest Strehl ratio. Since circumstellar signals are recovered for a range of |$n_{\rm pc, A}$| values, while the structure of residual speckle noise varies with |$n_{\rm pc, A}$|⁠, we median-combine the results obtained with |$n_{\rm pc, A}$| ranging from 1 to 8 and from 4 to 10 for the cropped and dither-combined cubes, respectively. Finally we combined into a single image, shown in Fig. 1, the results obtained with PCA-SADI on the cropped and dither-combined cubes in the inner and outer |$\sim$|0.15 arcsec radius, respectively.

3 RESULTS

3.1 Moving feature in a gap

Fig. 1 shows our reductions of all data sets including the new 2021 May 17 SPHERE/IFS data set, and JWST/NIRCam data presented in Christiaens et al. (2024). We re-detect PDS 70 b in all data sets, whereas PDS 70 c is best detected at |$\lambda >2~\mu\mathrm{m}$| due to the shape of its SED (e.g. Wang et al. 2020; Christiaens et al. 2024). The ‘PLF’ or ‘PDS 70 d’ is also re-detected in all observations at an average separation of |$\sim$|110 mas to the north-west of the star.

The source features an extended structure in observations using the RDI strategy. Given the point-like nature of planet b with the RDI strategy, the resolved emission may be genuine and trace authentic circumstellar signals in the direct vicinity of candidate d (mixed with some contamination from inner disc signal for the case of JWST). When performing only PCA-ADI on the 2022 February 28 star-hopping epoch (i.e. not using the reference star), we retrieved a more point-like morphology for d in YJH consistent with the earlier SPHERE epochs. However, this observation featured significantly less parallactic angle variation, exacerbating the self-subtraction of extended signals. Without a star-hopping epoch prior to 2021 it is unclear whether this is the true long-term morphology of the ‘PLF’ and its compact appearance was simply due to the geometric biases of ADI (Milli et al. 2012). We do, however, notice a dependence on the amount of parallactic angle variation (i.e. 2018 February 25, where it is slightly extended compared to shorter sequences).

In order to retrieve the astrometry of candidate d in the SPHERE data, we used the Negative Fake Companion (NEGFC) injection technique (Lagrange et al. 2010) as implemented in vip (described in appendix D of Christiaens et al. 2021). First, a Nelder–Mead simplex algorithm was used to obtain a guess on the negative planet flux and its location using |$\chi ^2_r = \sum \frac{(I- \mu)^2}{\sigma ^2}$|⁠, where I are the pixel intensities measured in a 1.5 |$\times$| FWHM aperture at the approximate location of candidate d, and |$\mu$| and |$\sigma$| the mean and standard deviation of pixel intensities in a 3|$\times$|FWHM-wide annulus at the separation of d. This estimate was provided to a Markov Chain Monte Carlo (MCMC) algorithm to sample the probability distribution of the (negative) flux, separation and PA, where |$\sigma$| was used for scaling the Gaussian log likelihood expression in the MCMC. To check for convergence we used a test based on the integrated autocorrelation time |$\tau$| as recommended in the emcee documentation, where the number of samples N must be significantly greater than |$\tau$|⁠. We used N/max(⁠|$\tau$|⁠) |$>$| 50 for each parameter leading to around 2500 iterations per data set. We then considered a 30 per cent burn-in on the MCMC chains. We accounted for the mean radial transmission of the apodized Lyot-stop coronagraph at the location of the injected negative companions (e.g. |$\sim$|30  per cent attenuation for a candidate at 110 mas as in the SPHERE manual). We used PCA-RDI for the star-hopping epochs and annular PCA-ADI otherwise. The SPHERE/IFS cubes were median combined in the spectral dimension to improve the SNR, and we considered a 500 |$\times$| 500 mas subframe centred on the star (for computational purposes) to build the PSF library. The relative astrometry obtained by the MCMC-NEGFC injection is given in Table 2. Our final reported astrometric errors are the quadratic sum of a Gaussian fit to the posterior distribution of separation and PA (Wertz et al. 2017), true North uncertainty (and 0|${_{.}^{\circ}}$|1 instrument offset uncertainty for the IFS), plate scale uncertainty, and centring uncertainty of 5 mas for coronagraphic sequences.

Table 2.

Relative astrometry of the candidate PDS 70 d around PDS 70 A obtained with the NEGFC technique described in Section 3.1.

Date	Separation (mas)	|$\sigma _{\rm sep}$| (mas)	PA (⁠|$^\circ$|⁠)	|$\sigma _{\rm PA}$| (⁠|$^\circ$|⁠)
2014-05-10	117.2	12.5	348.7	4.4
2016-06-02	114.3	17.1	334.6	11.7
2018-02-25	117.9	6.0	317.7	2.7
2019-03-06	120.4	8.3	315.7	2.9
2021-05-17	102.4	12.1	302.6	5.6
2021-08-21	110.7	12.7	305.7	6.5
2021-09-03	106.4	10.7	299.7	6.2
2022-02-28	109.3	6.5	299.4	3.9
2023-03-08	103.4	23.2	293.0	12.7

Notes. Final errors are a Gaussian fit to the distributions for each parameter, added quadratically with the centring uncertainty (5 mas for coronagraphic observations with SPHERE), plate scale error, and true North error.

As the NaCo data set contained residual signal at the same separation as d in the post-processed images, we opted to use a different technique than for SPHERE to better account for the biases of this noise and lack of pixels to estimate the noise. We used a Nelder–Mead simplex-NEGFC to identify the negative flux that best fits the curvature of the pixel intensities around d, by minimizing the absolute value of the determinant of the Hessian matrix in a 3 |$\times$| 3 pixel subframe (Quanz et al. 2015; Christiaens et al. 2024). We removed planet b from the pre-processed cube using the negative flux also found using this technique, and provided the subtracted cube to the speckle_noise_uncertainty function of vip. We used this function to inject fake companions at the same separation and flux as estimated for d, but at all azimuths in increments of 1|$^\circ$| excluding 10|$^\circ$| on either side of d, to estimate their retrieved astrometry. This approach provides a more reliable estimate of the uncertainties associated with residual speckle noise based on the distribution of differences between the injected flux and the retrieved flux for a planet at the same on-sky separation as d. Our final reported astrometric errors in Table 2 for this instrument are the quadratic sum of a Gaussian fit to the distribution of deviations on separation and PA from the injected values, true North uncertainty, and plate scale uncertainty.

3.2 Orbital predictions

We used Octofitter⁴ v6.0.0 (Thompson et al. 2023) to infer the posterior distribution of the orbital parameters for three simulated planets. We sample the posterior using stabilized variational non-reversible parallel tempering (Syed et al. 2021; Surjanovic et al. 2022), implemented in Pigeons.jl (Surjanovic et al. 2023). Our choice of package was primarily motivated by computation time required to reach convergence and the ability to fit multiple planet orbits simultaneously, as Octofitter is the fastest orbital fitting code publicly available. We used eccentricity (e), inclination (i), semimajor axis (a), argument of periastron (|$\omega$|), longitude of the ascending node (|$\Omega$|), and epoch of periastron, which is derived from the position angle at a reference epoch of MJD  = 58889, to define the orbit. For the b and c planets, we used literature astrometry summarized in Wang et al. (2021), with additional JWST/NIRISS, JWST/NIRCam, and MagAO-X measurements (Christiaens et al. 2024; Blakely et al. 2025; Close et al. 2025, respectively).

We fit for all four bodies (A, b, c, and d) simultaneously and include planet–star and planet–planet interactions. To explore the parameter space for d, we set unconstrained distributions on all priors except for separation and mass, which were set to 1–20 au and 0–10|$M_{\rm Jup}$|⁠, respectively. For the other planets, we implemented a prior on the mutual inclinations between b and the disc, c and the disc, and b and c, to be within |$\pm 10^{\circ }$|⁠. For the disc, we set i  = 128.3|$^{\circ }$| and |$\Omega = 156.7^{\circ }$| (Keppler et al. 2019). A Gaussian prior centred at 0.88|${\rm M}_{\odot }$| with a standard deviation of 0.09|${\rm M}_{\odot }$| was used for the stellar mass, using the measurement derived by Wang et al. (2021), including the same conservative 10  per cent uncertainty as was used in their orbital analysis. Finally, we included Gaussian priors on the parallax and proper motion of the star from Gaia DR3 (Gaia Collaboration 2023). A summary of all priors is provided in Table 3. We also imposed that the orbits of b, c, and d are non-crossing, and the semimajor axis of c could not exceed 40 au based upon modelling of the outer dust ring in Zhou et al. (2025). Fig. 2 provides the orbital fit results for the three protoplanets obtained for the maximum likelihood sample, along with 100 randomly drawn posterior samples. Fig. 3 shows the orbits retrieved for 1000 randomly drawn samples, and the corner plot is provided in Fig. D1.

Figure 2.

The maximum likelihood sample (solid) and 100 randomly drawn samples (transaprent) from the octofitter orbital solutions for three protoplanets PDS 70 b (blue), PDS 70 c (orange), and candidate d (green), without constraining to any resonances or enforcing coplanarity. Top: orbits as viewed on-sky, with our measured astrometry and associated error. Bottom: separation and PA against time.

Figure 3.

Overview of the PDS 70 system, including 1000 random samples from the orbital posteriors of planets b (blue curves), c (yellow), and candidate d (green). The priors and median posterior for each orbital element retrieved with octofitter are summarized in Table 3. A shift in declination of |$-0.02\,\,\mathrm{ arcsec}$| was applied to the continuum emission (colourscale) to enforce the c_smm blob to overlap with near-IR astrometry of the same epoch. The grey beam of 44 |$\times$| 37 mas with |$\mathrm{PA} = 51{_{.}^{\circ}} 5$| in the bottom-left is for the continuum emission originally presented in Benisty et al. (2021).

Table 3.

Octofitter orbital parameters and priors.

Element	Prior	Prior range	Median posterior
|$a_b$| (au)	Log Uniform	1–30	|$23.7^{+2.1}_{-2.0}$|
|$e_b$|	Uniform	0–0.99	0.12|$^{+0.07}_{-0.07}$|
|$i_b$| (⁠|$^{\circ }$|⁠)	Sine	|$-180$|–|$+180$|	|$128.1^{+3.7}_{-2.8}$|
|$\Omega _b$| (⁠|$^{\circ }$|⁠)	Uniform	0–360	171|$^{+6}_{-7}$|
|$\omega _b$| (⁠|$^{\circ }$|⁠)	Uniform	0–360	|$110^{+47}_{-46}$|
|$\theta _{\textrm {b}}$|	Uniform	0–360	142.17|$^{+0.03}_{-0.03}$|
|$M_{\textrm {b}}$| (⁠|$M_{\rm Jup}$|⁠)	Uniform	1–10	5.6|$^{+3.0}_{-3.0}$|
|$a_c$| (au)	Log Uniform	1–45	33.6|$^{+2.4}_{-2.6}$|
|$e_c$|	Uniform	0–0.99	0.05|$^{+0.05}_{-0.03}$|
|$i_c$| (⁠|$^{\circ }$|⁠)	Sine	|$-180$|–|$+180$|	131.7|$^{+3.2}_{-3.1}$|
|$\Omega _c$| (⁠|$^{\circ }$|⁠)	Uniform	0–360	155|$^{+5}_{-6}$|
|$\omega _c$| (⁠|$^{\circ }$|⁠)	Uniform	0–360	185|$^{+136}_{-156}$|
|$\theta _{\textrm {c}}$|	Uniform	0–360	277.3|$^{+0.1}_{-0.1}$|
|$M_{\textrm {c}}$| (⁠|$M_{\rm Jup}$|⁠)	Uniform	1–10	5.4|$^{+3.1}_{-3.1}$|
|$a_d$| (au)	Log Uniform	1–20	12.9|$^{+2.0}_{-2.1}$|
|$e_d$|	Uniform	0–0.99	0.16|$^{+0.16}_{-0.12}$|
|$i_d$| (⁠|$^{\circ }$|⁠)	Sine	|$-180$|–|$+180$|	146|$^{+12}_{-9}$|
|$\Omega _d$| (⁠|$^{\circ }$|⁠)	Uniform	0–360	177|$^{+145}_{-103}$|
|$\omega _d$| (⁠|$^{\circ }$|⁠)	Uniform	0–360	181 |$^{+110}_{-122}$|
|$\theta _{\textrm {d}}$|	Uniform	0–360	310.1|$^{+1.5}_{-1.6}$|
|$M_{\textrm {d}}$| (⁠|$M_{\rm Jup}$|⁠)	Uniform	0–10	5.2|$^{+3.3}_{-3.5}$|
|$M_{\textrm {*}}$| (⁠|$M_\odot$|⁠)	Normal	|$0.88 \pm 0.09$|	0.92|$^{+0.07}_{-0.07}$|
Parallax (mas)	Normal	|$8.897\pm 0.019$|	8.90|$^{+0.02}_{-0.02}$|

Notes. A Gaussian co-planarity prior of |$\pm 10^{\circ }$| was placed on the mutual inclination of b with the disc, c with the disc, and b with c. We imposed non-crossing orbits for all planets, and a constraint on c to not exceed 40 au. |$\theta$| is the position angle at a reference epoch of MJD |$=$| 58889.

For the known protoplanets b and c, we find similar, but notably different, parameters to those reported by Wang et al. (2021). We find a smaller stellar mass, along with a larger semimajor axis and lower eccentricity for planet b. To test if this can be explained by the presence of planet d, we compared our astrometry model with that of Wang et al. (2021) by testing a two-planet model. We found near-identical numerical results between orbitize! and octofitter for individual model evaluations, but were not able to replicate the same overall posterior distribution of the two-planet solution with either octofitter with pigeons or orbitize! with ptemcee. A future work could further investigate the source of this disagreement.

We were able to constrain the semimajor axis of PDS 70 d to 12.9|$^{+2.0}_{-2.1}$| au when including all four bodies in the fit without enforcing any orbital resonances. Inclinations approximately equal to that of the outer dust continuum ring (⁠|$\sim$|130|$^\circ$| in Casassus & Cárcamo 2022) feature the strongest likelihood within our prior range for b % c. The inclination and eccentricity of d, however, is not as constrained. Without the precision of VLTI/GRAVITY astrometry these orbital parameters will remain difficult to narrow down.

3.3 Spectrum

We opted to use the star-hopping strategy to retrieve the spectrum for the YJH and Ks bands as this approach results in the fewest geometric biases at short separation from the star (Wahhaj et al. 2021; Xie et al. 2022). For simplicity we consider only the 2018 February 25 epoch in K|$_1$|K|$_2$| as it was the longest sequence of its kind and was taken in the most favourable conditions.

For IFS, we used the following technique to retrieve the spectrum: (i) measure the SNR of candidate d in individual spectral channels (Fig. E1) to determine the 10 most favourable channels for astrometry, (ii) repeat the MCMC-NEGFC technique as in Section 3.1 after collapsing the 10 best channels to infer the separation and PA of d, (iii) run a Nelder–Mead simplex-NEGFC over all 39 spectral channels individually by fitting the curvature of the pixel intensities in a 5 |$\times$| 5 pixel subframe (as for NaCo, Section 3.1) to retrieve the contrast of d, with the separation and PA fixed to the value found in step (ii), and (iv) inject fake companions at the same separation and flux as estimated for d, but at all azimuths in increments of 1|$^\circ$|excluding 10|$^\circ$|on either side of b and d, and estimate their retrieved flux after PCA-RDI using Nelder–Mead simplex-NEGFC with the speckle_noise_uncertainty function of vip. The PCA-RDI images for each IFS channel after subtracting our best flux estimate for d is shown in Fig. E2. For IRDIS K|$_1$|⁠, K|$_2$|, and Ks bands we directly took the resulting flux from the MCMC-NEGFC technique in Section 3.1, and performed the speckle noise uncertainty procedure to retrieve errors on the flux.

For NaCo we used the flux and associated error recovered from the Nelder–Mead simplex-NEGFC and speckle noise uncertainty, as described in Section 3.2. We employed the same technique for SINFONI, however identified the best-fitting negative flux by minimizing the standard deviation of residual intensities in a 0.7|$\times$|FWHM aperture around d rather than fitting the curvature of the pixel intensities.

To convert from detector units to physical units, we took the contrast ratio of d and that of PDS 70 A measured in the non-coronographic PSF frames, and multiplied it by the flux-calibrated spectrum of PDS 70 A acquired with SpeX at the NASA Infrared Telescope Facility (IRTF). This choice of spectrum was motived by the need to include unresolved hot inner disc signal. As the SpeX data do not extend to 3.78 |$\mu\mathrm{m}$|⁠, we interpolated NASA/WISE W1 (3.35 |$\mu\mathrm{m}$|⁠) and W2 (4.6 |$\mu\mathrm{m}$|⁠) photometry in log space. For NIRCam we simply used the flux-calibrated F187N measurement of the star published in Christiaens et al. (2024).

There are two caveats with the above technique. First, d is at a radial separation that is partially attenuated by the SPHERE coronagraph by |$\sim$|30  per cent and this attenuation varies with coronagraph azimuth (and consequently, by epoch). We used the transmission profile for the N_ALC_YJH_S with the YJ prism as given in the SPHERE manual to scale the flux from d for each wavelength, which is not directly equivalent to the radial transmission expected for the YJH prism. Secondly, PDS 70 A was placed behind the coronagraphic mask for the FLUX cubes during the 2022 February 28 star-hopping epoch, meaning there was no direct observation of this star on this night. We instead used FLUX cubes taken in similar observing conditions and with the YJH prism from the 2021 May 17 observation for IFS and the 2021 September 3 observation for IRDIS, ensuring to correctly scale the DIT. Any variability of the star at these wavelengths between these observations would be folded into the spectrum.

Fig. 4 shows the spectrum and associated uncertainty. The shape is consistent with the SED of the star and unresolved inner disc spectrum at short wavelengths for the RDI strategy; however the uncertainty due to speckle noise is significant at the separation of the source. The NIRCam and IRDIS measurements follow the same trend as the IFS, whereas the longer wavelength SINFONI and NaCo measurements are significantly brighter.

Figure 4.

Extracted spectrum of the PDS 70 d candidate for all instruments presented in this work. For the 0.95–1.65 |$\mu\mathrm{m}$| (blue) and 2.18 |$\mu\mathrm{m}$| (orange) measurements we used the 2022 February 28 star-hopping epoch, while the 2.11 and 2.23 |$\mu\mathrm{m}$| (orange) measurements are from 2018 February 25. The resampled stellar spectrum (gray) we used to convert to flux units is from SpeX at the IRFT and was scaled by 3|$\times 10^{-4}$| to match the intensity of candidate d at recent epochs. The horizontal error bars correspond to the effective filter width for each instrument. The spectral slope of the candidate is overall consistent with that of the star, except for the SINFONI and NaCo points for which possible explanations are discussed in text.

Several mechanisms can produce the high flux we measured at 2.37 and 3.78 |$\mu\mathrm{m}$|⁠. PDS 70 A is known to be highly variable (Perotti et al. 2023; Gaidos et al. 2024; Jang et al. 2024) likely due to occulting material close to the star. This is also indicated by the WISE W1 (3.35 |$\mu\mathrm{m}$|⁠) time-series observations, which is the closest WISE band to our L’ measurement. If the star and unresolved inner disc were dimmer during the SINFONI and NaCo epochs, we would overestimate the 2.37 and 3.78 |$\mu\mathrm{m}$| flux. For this reason, we corrected all of our flux measurements for the stellar variability measured in WISE W1 by assuming a consistent correlation for all wavelengths. The largest change applied to the SINFONI epoch, where PDS 70 A dimmed by |$\sim$|40  per cent at 3.35 |$\mu\mathrm{m}$| (see fig. 7 of Gaidos et al. 2024). This may none the less not be a sufficient correction given the sparse sampling of stellar variability and the assumed linear relation between measurements, which may underpredict the flux drop at the SINFONI epoch if already undergoing the dipper event. Another source of bias is the inclusion of an unknown amount of contamination from the inner disc in the PSF. This is most relevant for NaCo where the FWHM during the sequence was 108 mas, although may also be relevant to the SINFONI data where the achieved Strehl ratio is likely lower than in the SPHERE data.

4 DISCUSSION

4.1 Evidence for a third protoplanet

The source is detected across nine epochs over 9 yr. This lends support to it being true signal rather than a post-processing artefact, and can rule out the hypothesis of the signal tracing a static inner disc or residual speckle. A speckle would be transient across observations and not be consistent through the instruments included in this work. Candidate d is also detected in individual spectral channels (Appendix  E). Additionally, the source is not in the predicted location of the Lagrange points L|$_4$| or L|$_5$| for either b or c (Balsalobre-Ruza et al. 2023). It is therefore not expected to be tracing a dust trap related to these protoplanets. Additionally, mayonnaise identifies it as a point source rather than a component of the disc (Appendix  B), although there is another feature at the same separation pointed out as PLF 2 in Pairet et al. (2021). Note that we also detected d in the IRDIS 2018 Feb 25 data with the andromeda (Cantalloube et al. 2015) post-processing algorithm, which filters out extended structures to find point-like sources (Appendix  B). While not terribly revealing, a third body easily fits within the observed astrometric Reduced Unit Weight Error (RUWE) reported from Gaia (Appendix  C).

Orbital fitting in Section 3.2 shows that d can be described as a bound object with Keplerian frequency approximately in the plane of the disc – as are b and c. Interestingly, the favoured solutions allude to a Laplace resonance with the two known protoplanets. Considering the uncertainties on the period of candidate d and on the stellar mass, the exact resonance is not clear. The median posterior from octofitter gives a period ratio of d and b as 2.46|$^{+0.82}_{-0.53}$|⁠, also including the possibility of the 3:1 or 5:2 resonance, and of c and d it is 4.16 |$^{+1.32}_{-0.93}$|⁠. Using the dynamically stable result for b in Wang et al. (2021) gives a period ratio of d and b of 2.05|$^{+0.74}_{-0.48}$|⁠. Given the mass of the other bodies we expect such an MMR in order for a chain of three giant planets to be stable over long time periods (Myr). As the age of PDS 70 A is 5.4|$\pm$|1 Myr (Müller et al. 2018) any protoplanets would need to be in strong, preferably first-order, resonance to avoid ejection. Such a scenario would be similar to the directly imaged giant planets HR 8799 b, c, and d (Goździewski & Migaszewski 2020), the Gliese-876 b, c, and d giant planets (Nelson et al. 2016; Cimerman et al. 2018), and the Galilean moons, which all share a Laplace resonance.

An explanation for the SINFONI flux could be the first CO overtone v  = 0 to v  = 2, which occurs at 2.38 |$\mu\mathrm{m}$|⁠. Such excitation would require warm gas in the vicinity of a planet, heated by stellar radiation (Oberg et al. 2020), by radiation from the planet, and by accretion. Similar predictions have been made for |$^{12}$|CO rotational–vibrational lines in Oberg et al. (2023). Alternatively, it may probe thermal excess towards the end of the K band from a circumplanetary disc (CPD). Interestingly, the flux estimated in the NaCo data lends support to the CPD hypothesis. This would also be consistent with b and c where circumplanetary contribution is required to match the SED of the planets at these wavelengths (Christiaens et al. 2019b; Blakely et al. 2025; Choksi & Chiang 2025).

The extended nature of d in the RDI observations is suggestive of a spiral wake or circumplanetary envelope, if d is a protoplanet. A protoplanet would dynamically interact with the inner disc and release spiral density waves (Ogilvie & Lubow 2002). In this scenario, the signal seen at short wavelengths could therefore be starlight scattering off small grains in the outer wake of the protoplanet.

4.2 Evidence against a third protoplanet

If indeed d is a protoplanet at a deprojected separation of 12.2 au, it is in the same annular gap as the other protoplanets. However, the extent of the inner disc is not well constrained and needs to be considered. Recent ALMA 855 |$\mu\mathrm{m}$| observations in Benisty et al. (2021) show a mostly resolved inner disc, and at higher resolutions its distribution is irregular. The ALMA data show the presence of mm-emitting grains close to the separation of candidate d. However the orbital distance of the mm-emitting grains shrinks with a decreasing robust parameter (robust  = 0.5 in our Fig. 3), and an inner disc is not detected in the 3 mm continuum (Doi et al. 2024). Benisty et al. (2021) modelled the 855 |$\mu\mathrm{m}$| signal using a Gaussian with a width of 59 mas. By taking 2|$\sigma$| (twice the width of the Gaussian) we get an inner disc size of about 120 mas, which is 60 mas each side of the star (7 au projected). Candidate d is exterior (⁠|$\sim$|110 mas) to such an inner disc assuming a uniform distribution of grains.

By paying special attention to uncertainties and the JvM correction (Jorsater & van Moorsel 1995) of the same 855 |$\mu\mathrm{m}$| data, Casassus & Cárcamo (2022) further unveiled the finer structure of the inner disc. It is highly irregular and the morphology varies significantly over three epochs (2016, 2017, and 2019). Most notably, there is a detached clump of emission towards the North at 0.08 arcsec, or 4.5 au projected, separation from the ring centre. Such structures seem to change between epochs, indicating they are either variable or artefacts that require further study. We stress that such emission is interior to d observed in the near-IR, although the presence of variable, ambient material in proximity complicates any claim that d could be a protoplanet and detached from the inner disc. If it is physical, what drives this irregular structure is unclear, but could be related to the gravitational perturbation from a nearby protoplanet.

Many of the orbital solutions in Fig. 3 pass close to the 855 |$\mu\mathrm{m}$| inner disc continuum with a beam of 44 |$\times$| 37 mas. A protoplanet of similar mass to PDS 70 b and c should have carved out a dust gap in its orbital path. Whether the irregular and clumpy distribution of the inner disc can be attributed to dynamical interaction with a protoplanet at |$\sim 13$| au, or simply to b alone, is difficult to disentangle. New hydrodynamical simulations with the addition of a third protoplanet at the separation of d could indicate whether the inner disc is stable over several Myr with such a configuration of protoplanets.

4.3 Implications for future observations

Based on our fits for a Keplerian protoplanet we expect planet d to orbit clockwise and reach the semiminor axis of the disc within a few years. For this reason, its polarized intensity will decrease due to the scattering angle dependence (forward versus backward scattering), whereas the total intensity will increase. While this is favourable for future studies, the predicted orbits indicate that d will pass further behind a coronagraphic mask due to the viewing angle on-sky. This caveat may motivate non-coronagraphic follow-up campaigns, particularly with star-hopping, to reach deeper contrast at short separation. Such an observation would need to be taken in more stable observing conditions than previously due to the presence of bright speckles at close separations.

Follow-up observations with VLTI/GRAVITY are an obvious next step to significantly constrain the orbit of d, as has already been demonstrated with two epochs on both b and c (Wang et al. 2021). The detection of d in K bands with SPHERE/IRDIS suggests that it will be detectable in the operating wavelengths of GRAVITY (2.0–2.4 |$\mu\mathrm{m}$|⁠) and such an observation would provide a higher resolution spectrum. The |$\sim$|50 mas field of view of GRAVITY is wide enough to capture all orbital solutions found with existing HCI in this work for at least the next few years. However the presence of bright extended signals, such as those present in the star-hopping data in the vicinity of d, could complicate an observation with GRAVITY.

Thermal emission, or lack thereof, in ELT/METIS observations may be required to fully distinguish between the hypothesis of d being a disc signal or a third protoplanet, as thermal contribution is expected to come either directly from the surface of the planet or from the CPD (Szulágyi et al. 2019; Choksi & Chiang 2025). The higher angular resolution of the ELT at these wavelengths would disentangle inner disc signal from any real planet signal and enable more robust spectral characterization.

5 CONCLUSIONS

We present 9 yr of high-contrast imaging data of PDS 70 with VLT/SINFONI, VLT/NaCo, VLT/SPHERE, and JWST/NIRCam to obtain estimates on the orbit and spectrum of a third protoplanet candidate. Our analysis of 10 data sets over eight epochs plus the published JWST/NIRCam observations suggest the following:

• The candidate protoplanet PDS 70 d is detected at 0.95–3.8 |$\mu\mathrm{m}$|⁠, irrespective of observing strategy (ADI, SDI, or RDI) or post-processing technique, in clockwise motion consistent with Mesa et al. (2019) and Christiaens et al. (2024).

• Orbital solutions suggest a near 4:2:1 Laplace resonance with the existing protoplanets in the plane of the disc, although the eccentricity is uncertain. The semimajor axis of d is 12.9|$^{+2.0}_{-2.1}$| au when including all four bodies in the fit (A, b, c, d).

• The YJH spectrum of d is consistent with the SED of the star and unresolved inner disc, suggestive of starlight scattering off small dust grains and is also bluer than planet b.

• A thermal excess is detected at 2.37 and 3.78 |$\mu\mathrm{m}$|⁠. Its origin is complicated by stellar variability and contamination from the inner disc, but such an excess can come from a CPD.

• The source may not be ‘point-like’ as previously reported, and instead appears extended in RDI data due to the absence of ADI self-subtraction. To explain such morphology, the feature could be tracing an outer spiral wake from a third protoplanet with a dusty envelope, and could explain the lack of accretion signatures.

• The extent of the inner disc around PDS 70 A is not well constrained and, in some cases, is close to the predicted orbits of PDS 70 d.

ACKNOWLEDGEMENTS

We thank Myriam Benisty for useful discussions and for kindly supplying the ALMA data. We also thank Yifan Zhou for discussions regarding the HST observations. IH acknowledges a Research Training Program scholarship from the Australian government. VC and SJ thank the Belgian F.R.S.-FNRS. VC thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000142531. GDM acknowledges the support from the European Research Council (ERC) under the Horizon 2020 Framework Program via the ERC Advanced Grant ‘ORIGINS’ (PI: Henning), Nr. 832428, and via the research and innovation programme ‘PROTOPLANETS’, grant agreement Nr. 101002188 (PI: Benisty). We made use of the Multi-modal Australian ScienceS Imaging and Visualization Environment (MASSIVE) and data from the European Space Agency (ESA) mission Gaia. This work was performed on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government, and from the Victorian Higher Education State Investment Fund (VHESIF) provided by the Victorian Government. We are grateful for Australian Research Council Discovery Project funding via DP220103767 and DP240103290.

The scientific results presented in this article made use of Python and the numpy, scipy, astropy, and matplotlib packages.

DATA AVAILABILITY

Data are available as FITS files through the Max Planck Digital Library’s Keeper research repository via the following link: https://keeper.mpdl.mpg.de/d/2633616c51b64433962d/

Footnotes

1

https://github.com/VChristiaens/vcal_sphere

2

https://github.com/bpairet/mayo_hci

3

https://github.com/PynPoint/PynPoint

4

https://sefffal.github.io/Octofitter.jl

APPENDIX A: GEMINI/NICI L’-BAND ANGULAR DIFFERENTIAL IMAGING OBSERVATIONS

We reprocessed non-coronagraphic Gemini/NICI observations from the night of 2012 March 31 to extend our time baseline to almost 10 yr. These L’-data, in pupil stabilized mode, were published in Hashimoto et al. (2012). Details on the pre-processing are in Keppler et al. (2018). The final cube consisted of 107 frames (comprised of 25 co-adds of 0.76s) with 99.4° of parallactic angle variation after bad frame rejection. We used the same approach to post-processing as with NaCo, using annular PCA-ADI in a 1.2 arcsec wide annular centred on the star and a rotation threshold of 0.5|$\times$|FWHM at the approximate separation of candidate d and explored |$n_{\rm pc}$|  = 1–50 with |$n_{\rm pc}$|  = 10 shown in Fig. A1. A PLF is identified at |$\sim$|336.4° |$\pm$| 10.87° with a radial separation of 146.99 |$\pm$| 19.91 mas. We suspect this feature to trace a residual speckle rather than thermal emission from d, due to its comparatively wider separation that overlaps with the predicted orbital solutions for b. The feature is also of low significance given other brighter signals at the same separation.

Figure A1.

Our post-processed image of the Gemini/NICI data from 2012 March 31 displayed the same as all other epochs in Fig. 1. Several signals are retrieved at the same approximate separation as PDS 70 d.

APPENDIX B: ANDROMEDA AND MAYONNAISE REDUCTION OF SPHERE DATA

The SNR maps obtained with andromeda for the 2018 February 25 SPHERE data are shown in Fig. B1. The two known protoplanets are detected by each subinstrument, whereas signal from candidate d is significantly more robust in the IFS wavelength range. The successful detection indicates that d traces a rotating point source. We also tested the mayonnaise algorithm (Pairet et al. 2021) on IRDIS K|$_1$| data to disentangle point sources from extended signals. Fig. B1 (right) shows the results of this test. The algorithm does identify d as well as b and c as non-disc features, but also does the same for emission to the east of the star (PLF 2, Pairet et al. 2021).

Figure B1.

Left and Middle: IFS and IRDIS (K|$_1$|K|$_2$|) SNR maps obtained with the andromeda algorithm. Right: mayonnaise (Pairet et al. 2021) identifies PDS 70 d as a point source and not a disc feature, but it does the same for the emission in the east of the star.

APPENDIX C: RUWE PREDICTIONS

As the candidate is closer to the star than the existing planets, we tested whether a third protoplanet is compatible with the re-normalized unit weight error (RUWE) of 1.35 from Gaia DR3 (Gaia Collaboration 2023). This value is calculated as a |$\chi ^2$| error by comparing the position of the photocentre to a single star solution (Penoyre et al. 2020) and re-normalized for a star’s brightness and colour (Lindegren 2018). An RUWE value of near 1 indicates the photocentre’s position can be explained by a single star track. Values greater than |$\sim$|1.25 could be indicative of a companion (Penoyre, Belokurov & Evans 2022).

In order to understand the effect of a planetary-mass companion on RUWE, we simulated the value of RUWE of PDS 70 A as a function of planet mass and period, assuming a single-planet model and a distance of 113 pc. Inclination was fixed at 130|$^{\circ }$| and we ignore the presence of the disc, assuming RUWE is solely influenced by a companion. This is presented in Fig. C1. This result demonstrates that the existing protoplanets, even at their highest mass estimates, have effectively no impact on the RUWE due to their separation from the star. Adding another protoplanet at 13 au does not significantly affect the predicted RUWE. The distance to the source and other factors such as variability (Gaidos et al. 2024) and obstructing disc material are likely dominating the RUWE budget more than any of the proposed or confirmed planets.

Figure C1.

RUWE of PDS 70 A as a function of planet mass and period. This demonstrates that the detected planets, or PDS 70 d, cannot alone account for the measured RUWE of 1.35.

APPENDIX D: OCTOFITTER CORNER PLOT

In Fig. D1 we show the marginal posterior distributions of the eccentricity, mutual inclination with the disc, and orbital period of PDS 70 b, PDS 70 c, and PDS 70 d, from our model described in Section 3.2.

Figure D1.

Comparison between the posterior samples of the eccentricity, mutual inclination with the disc, assuming |$i = 128{_{.}^{\circ}} 3$| and |$\Omega = 156{_{.}^{\circ}} 7$| for the disc (Keppler et al. 2019), and orbital period of the protoplanets b and c and candidate d, from our Octofitter results, described in Section 3.2.

APPENDIX E: INDIVIDUAL IFS CHANNELS

Fig. E1 shows all 39 SPHERE/IFS channels from the 2022 February 28 epoch after subtraction of the reference star. Candidate d is detected in channels that do not suffer from telluric absorption, adding support to it being a genuine feature rather than an observational artefact. Fig. E2 shows the same cube with the best estimate for the flux of d subtracted.

Figure E1.

IPCA-RDI for the 2022 February 28 epoch in Fig. 1 truncated to seven principal components and obtained with three iterations, for each individual spectral channel. For all frames, the PSF of the reference star UCAC2 14412811 observed during the sequence has been subtracted. Candidate d is visible in all channels that are not contaminated by telluric absorption, whereas PDS 70 b is conspicuous at longer wavelengths (i.e. H band). The grey patch represents the FWHM of the non-coronagraphic PSF and is unique for each channel. No numerical mask has been used to represent the location of the coronagraphic mask. Pixel intensities are on a linear scale with an intensity cut from the 0–99.5 percentiles of the entire cube. Figure made using hciplot (https://github.com/carlos-gg/hciplot).

Figure E2.

Same as Fig. E1, but with our best estimate for the flux of PDS 70 d from Section 3.3 for each channel subtracted. The flux is generally well subtracted with the exception of some extended emission in the vicinity of d at short wavelengths.