Accretion Properties of PDS 70b with MUSE^*

Observations of PDS 70 were made with the optical integral field spectrograph, the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010), at the Very Large Telescope (VLT) on 2018 June 20 (UT) under the clear sky condition with optical seeing of 070–081, as part of the commissioning runs of the Narrow-field Mode. The adaptive optics (AO) system for MUSE consists of the four laser guide star (LGS) facility and the deformable secondary mirror (DSM). The LGS wave front sensors measure the turbulence in the atmosphere, and the DSM corrects a wave front to provide near diffraction limited images. The data used in this work (a field of view of 742 × 743, a spatial sampling of 25 × 25 mas², a spectral range of 4650–9300 Å, a resolving power of 2484 ± 5 at 6500.0 Å) are publicly available on the ESO archive.¹³ Six individual raw frames with an exposure time of 300 s obtained in the observations were calibrated with the MUSE pipeline v2.6 (Weilbacher et al. 2012; see Haffert et al. 2019 for more details about data). The absolute flux was also calibrated in the pipeline where the flux calibration curve is derived from both an atmospheric extinction curve at Cerro Paranal and a spectro-photometric standard star stored as MUSE master calibrations. The apparent magnitude of PDS 70 in calibrated MUSE data at 6905–8440 Å corresponding to the Sloan $i^{\prime}$ band filter is measured to be 12.0 ± 0.1 mag, which is fainter than the literature value of i' = 11.1 ± 0.1 mag (Henden et al. 2015). In the following post-processing, we used an image size of 40 × 40 spatial pixels (corresponding to 1'' × 1'') around the central star with a high signal-to-noise ratio (S/N). In changing the size, we found that this size maximized the S/N of PDS 70b. Furthermore, since there are no signals in ∼5780 to ∼6050 Å to avoid contamination by sodium light from LGS, we removed them from the data.

After the basic calibration with the MUSE pipeline via a scientific workflow of EsoReflex (Freudling et al. 2013), we extracted six data cubes and inspected them carefully. We found that five of the data cubes have a relatively strong stripe pattern in the vicinity of the central star (Figure A1 in Appendix), induced by a wavelength shift due to possible calibration errors. Since we do not see a wavelength shift in one data cube, we conducted cross-correlation analyses to correct the wavelength shift. The FWHM of the point-spread function (PSF) in the final image was 66 mas at 6562.8 Å.

To subtract the stellar spectra of PDS 70 at each spatial pixel, we followed the methods used in Haffert et al. (2019). An essential point is to construct reference spectra from the most correlated components obtained from principal component analysis (PCA; see the basis of PCA given by Jee et al. 2007). Before applying PCA, we normalized each spectrum by the median values of each spectrum (purple line in Figure A3 in the Appendix). We then calibrated a fake continuum pattern, as shown in Figure A2 in the Appendix, where we show that the stellar spectra are different from the halo spectra at 300 mas from the stellar position. This is because the degree of central concentration of the stellar flux is higher at longer wavelengths due to the better AO performance at longer wavelengths. To calibrate this fake continuum, we generated median spectra for 40 × 40 spatial pixels in six data cubes (blue line in Figure A3). We chose the median to avoid residual bad/hot pixels. After dividing each spectrum by the median spectrum (red line in Figure A3), we smoothed each divided spectrum by a Gaussian function with a kernel of 230 spectral channels (black line in Figure A3). Finally, we divided each normalized spectrum (purple line in Figure A3) by the smoothed spectrum (black line in Figure A3) to correct the fake continuum (yellow line in Figure A3).

We applied PCA to 9600 spectra (six frames with 40 × 40 spatial pixels) and found the optimum number of components to subtract to be 1 by maximizing the S/N of PDS 70b. After subtracting reference spectra for each spatial pixel, we recovered the unit (10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹) in the image by performing inverse processes of the above normalization. Two data cubes out of six were removed from the final image since they degraded the data quality, i.e., gave a worse S/N of PDS 70b. The total exposure time after removing the two cubes was 1200 s. The final images are shown in Figure 1, in which the Hα and Hβ images are constructed with three (656.096, 656.221, and 656.346 nm) and one (486.096 nm) wavelength channels, respectively. Astrometry of the two planets is summarized in Table 1.

Zoom In Zoom Out Reset image size

The Hα fluxes were measured with a box aperture of 3 × 3 pixels centered on PDS 70b and c in the combined image of three channels. The errors of Hα and Hβ per spectral channel were estimated from the standard deviation of the flux within the box aperture (3 × 3 pixels) at the planetary positions in a continuum spectra, i.e., ∼6570 to ∼6700 Å and ∼4800 to ∼4930 Å for Hα and Hβ, respectively. The S/Ns of Hα were 27 and 10 for PDS 70b and c, respectively (Table 1), while the peak S/Ns along the spectral direction were 20 and 9 for PDS 70b and c, respectively (Figure 1).

To recover the missing flux, we applied two flux calibrations. The first was an aperture correction (Howell 1990) in which we recovered the missing flux due to the small aperture. The other was a flux correction in which we estimated the lost flux during the PSF subtraction. In the aperture correction, the correction factors (between the box aperture of 3 × 3 pixels and a circular aperture of 70 pixels in radius) derived by the photometry of the primary star were 6.65 and 16.9 for Hα and Hβ, respectively. Since the AO performance is better at longer wavelengths, the correction factor of Hα is smaller than that of Hβ. We then inspected the extent to which the point-source flux decreased during the PSF subtraction. Positive fake sources were injected in the Hα and Hβ channels, and the PSF subtraction process was rerun. The fake sources were injected at the separations of PDS 70b and c, and at position angles of 0°–345° with 15° steps, except the planetary positions, and with fluxes corresponding to those of PDS 70b and c without flux calibrations. Regarding the Hβ flux, we used a 3σ value as the input flux. The results show that the flux decreased to 90%, 98%, 90%, and 95% for PDS 70b at Hα, PDS 70c at Hα, PDS 70b at Hβ, and PDS 70c at Hβ, respectively. The corrected line fluxes and their errors are shown in Table 1. The linewidths were measured by fitting a Gaussian profile to the Hα spectra with the fit function in GNUPLOT.

The simultaneous observations of Hα and Hβ with MUSE provide reliable measurements of the line flux ratio without any concerns regarding the temporal variability of each flux. Signals of Hβ emissions are not seen in both PDS 70b and c, and we only obtained the upper limit of its flux. However, the upper limit of Hβ emissions allows us to constrain the physical parameters related to the planetary accretion. Following the method in Aoyama & Ikoma (2019), we first estimated the free-fall velocity v₀ of the gas and the pre-shock number density of hydrogen nuclei n₀ using the accretion shock model in Aoyama et al. (2018). Here, we assume that the accretion speed is the same as the free-fall velocity. The observational constraints for these quantities are given by the Hα linewidths at 10% and 50% maxima (Table 1). As explained in Aoyama & Ikoma (2019), the Hα linewidth increases with n₀ and v₀. This is because an increase of n₀ enhances the absorption of Hα emissions near the line center by shock-heated gas, while a high value of v₀ means that a flow with higher velocity passes. Thus the Hα linewidth is a function of n₀ and v₀ (see Figure 5 in Aoyama & Ikoma 2019). In Figure 3, the intersections of the Hα 10% and 50% linewidths indicate that (v₀, n₀) is roughly (144 km s⁻¹, 3.8 × 10¹³ cm⁻³) and (139 km s⁻¹, 3.2 × 10¹³ cm⁻³) for PDS 70b and c, respectively.

Zoom In Zoom Out Reset image size

The theoretical model of Aoyama & Ikoma (2019) suggests that the flux ratio of Hβ to Hα (${F}_{{\rm{H}}\beta }/{F}_{{\rm{H}}\alpha }$) increases with n₀ when there is no foreground extinction, e.g., interstellar extinction (right panel in Figure 3). This is because the higher density of n₀ ≳ 10¹³ cm⁻³ is large enough to saturate the Hα flux due to the effect of absorption in the post-shock region. Meanwhile, since the emission source for Hβ is more excited than that for Hα, higher values of v₀ or n₀ are necessary to put Hβ in such a saturated state (see Aoyama et al. 2018 for more details). Hence the value of F_Hβ/F_Hα is close to unity for n₀ ≳ 10¹³ cm⁻³. In contrast, a Hα absorption is negligible with the lower density of n₀ ≲ 10¹² cm⁻³, resulting in a smaller value of F_Hβ/F_Hα ≲ 0.5.

With values of n₀ and v₀ as derived above for PDS 70b and c, the theoretical value of the flux ratio is close to unity (right panel in Figure 3). However, our observed values of ${F}_{{\rm{H}}\beta }$/F_Hα at a 3σ upper limit for PDS 70b and c are <0.28 and <0.52, respectively, which are smaller than the theoretically expected values. These observed values are also similar with accreting T Tauri stars and brown dwarfs with the values of 0.1–0.5 (see Table 7 in Herczeg & Hillenbrand 2008). This underestimating of the flux ratio could be due to the extinction of circumplanetary materials. According to ALMA observations and numerical simulations by Keppler et al. (2019), gas is likely to exist in the dust gap of the circumstellar disk where PDS 70b and c are located. We conjecture that small dust grains (sub-micron size) are well coupled with gas in the gap, contributing to the extinction. Based on this speculation, we use the extinction law A_λ ∝ λ^−1.75 (Draine 1989) to estimate the extinction. Considering the theoretical prediction that the value of F_Hβ/F_Hα without extinction is unity, we derived A_Hα > 2.0 and >1.1 mag for PDS 70b and c, respectively. We note that the extinction law is valid in the range of 0.7 ≲ λ ≲ 5 μm, and the slope of extinction will be gentle at Hα and Hβ wavelengths (see Figure 1 in Draine 1989). The gentle slope will make the lower limit of A_Hα larger than the above value.

We estimate the filling factor ${f}_{{\rm{f}}}$ using the following equation for the Hα luminosity (Aoyama & Ikoma 2019):

$$\begin{eqnarray}&&{L}_{{{\rm{H}}}_{\alpha }}=4\pi {R}_{{\rm{p}}}^{2}{f}_{{\rm{f}}}{I}_{{{\rm{H}}}_{\alpha }}{10}^{-0.4{A}_{{{\rm{H}}}_{\alpha }}},\end{eqnarray} \tag{ 1 }$$

or

$$\begin{eqnarray}&&0.4{A}_{{{\rm{H}}}_{\alpha }}-{\mathrm{log}}_{10}{f}_{{\rm{f}}}=-{\mathrm{log}}_{10}{L}_{{{\rm{H}}}_{\alpha }}+{\mathrm{log}}_{10}(4\pi {R}_{{\rm{p}}}^{2}{I}_{{{\rm{H}}}_{\alpha }}),\end{eqnarray} \tag{ 2 }$$

where R_p and ${I}_{{{\rm{H}}}_{\alpha }}$ are the planetary radius and the Hα energy flux per unit area, respectively. The Hα luminosities are calculated from the line fluxes listed in Table 1. We assume that the planetary radius R_p is 2R_Jup, where R_Jup is the Jovian radius (Spiegel & Burrows 2012). ${I}_{{{\rm{H}}}_{\alpha }}$ is a function of n₀ and v₀, and its dependence is investigated in Aoyama et al. (2018). Using these results, we plot the lines of $0.4{A}_{{{\rm{H}}}_{\alpha }}-{\mathrm{log}}_{10}{f}_{{\rm{f}}}$ with different values in the v₀–n₀ plane (left and middle panels in Figure 3). The lower limits for the extinction A_Hα constrain the possible range of the filling factor: A_Hα > 2.0 and >1.1 mag give ${f}_{{\rm{f}}}\gtrsim 0.01$ and ≳0.003 for PDS 70b and c, respectively.

We can also derive a planetary dynamical mass M_P and an accretion rate ${\dot{M}}_{{\rm{P}}}$ (Aoyama & Ikoma 2019), written as

$$\begin{eqnarray}&&{M}_{{\rm{P}}}=\displaystyle \frac{{R}_{{\rm{P}}}{v}_{0}^{2}}{2G}\,\mathrm{and}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{P}}}=\mu {n}_{0}{v}_{0}{f}_{{\rm{f}}}(4\pi {R}_{{\rm{P}}}^{2}),\end{eqnarray} \tag{ 4 }$$

respectively, where G is the gravitational constant and μ is the mean mass per hydrogen nucleus. Using μ = 2.3 × 10⁻²⁴ g and R_P = 2 R_Jup, we estimate M_P of 12 ± 3 M_Jup and ${\dot{M}}_{{\rm{P}}}\,\gtrsim 5\times {10}^{-7}$ M_Jup yr⁻¹ for PDS 70b, and M_P of 11 ± 5 M_Jup and ${\dot{M}}_{{\rm{P}}}\gtrsim 1\times {10}^{-7}\,{M}_{\mathrm{Jup}}\,{\mathrm{yr}}^{-1}$ for PDS 70c.

The planetary mass is commonly estimated by comparing photometric and spectroscopic observations with planetary evolution and atmospheric models (e.g., Keppler et al. 2018; Müller et al. 2018; Christiaens et al. 2019a; Haffert et al. 2019; Mesa et al. 2019). On the other hand, our method based on the accretion shock model can estimate a dynamic planetary mass (Aoyama & Ikoma 2019). In Section 4.3, we estimated the masses of PDS 70b and c to be ∼12 and ∼11 M_Jup, respectively. As shown in Equation (3), these values depend on the planetary radii, which cannot be constrained only by MUSE observations. We consider the possible range of the planetary radii suggested by the planetary evolution model, which predicts the radii to be ∼1–2 R_Jup with an age of ≲10 Myr (Spiegel & Burrows 2012). Therefore, the planetary mass derived from the radii of 2R_Jup will give the upper limit. We also note that the observed linewidths will be the upper limits, because the linewidths of the two planets are similar to MUSE's spectral resolution, ∼120 km s⁻¹. Our estimated masses are consistent with previous near-infrared photometric and spectroscopic estimations: the mass of PDS 70b is estimated to be 5–9 and 2–17 M_Jup from photometric and spectroscopic observations, respectively (Keppler et al. 2018; Müller et al. 2018), while PDS 70c's mass is 4–12 M_Jup from photometric observations (Haffert et al. 2019; Mesa et al. 2019). Our mass estimation method is independent of the previous methods based on the planetary evolution and atmospheric models only, and therefore will provide a powerful tool to calibrate these models.

The extinction of the primary star was found to be negligible (Müller et al. 2018). However, our observations suggest that PDS 70b and c have A_Hα > 2.0 and >1.1 mag, respectively. We hypothesize that small (sub-micron size) grains coupled with the gas in the dust gap in the circumstellar disk are the cause of the planetary extinction. Such dust, if it exists, will not be clearly visible in existing observations because of their faintness. Here, we estimate the contribution of the unseen dust to the extinction. Keppler et al. (2019) conducted hydrodynamic simulations to reproduce the observed ¹²CO (2–1) integrated intensity map and found that the azimuthally averaged gas surface density at the location of PDS 70b is ∼0.01 g cm⁻², which can be translated to a molecular hydrogen surface density of ∼3 × 10²¹ cm⁻². By using the relationship between the interstellar dust extinction and the hydrogen column density (N_H/A_V = 1.9 × 10²¹ cm⁻² mag⁻¹; Bohlin et al. 1978), we derived A_V ∼ 3.3 mag, which corresponds to A_Hα ∼ 2.4 mag with A_λ ∝ λ^−1.75 (Draine 1989), where we assume that hydrogen atoms are in the form of hydrogen molecules. The derived vertical extinction at the gap region is comparable to the planetary extinction estimated by the accretion shock model, which supports our idea that the origin of extinction is small, unseen dust grains in the gap.

We now consider the influence of the extinction on the estimation of the mass accretion rate. First, we briefly explain the relationship between the extinction and the mass accretion rate. As shown in Equation (4), the mass accretion rate is a function of ${f}_{{\rm{f}}}$, n₀, and v₀. The values of n₀ and v₀ can be estimated from the Hα linewidths (left and middle panels in Figure 3), and hence we only need to estimate the value of ${f}_{{\rm{f}}}$. This value can be estimated from Equation (1), with the free parameter A_Hα. The value of A_Hα was independently estimated from the line ratio (right panel in Figure 3). For PDS 70b, we derived an Hα extinction of A_Hα > 2.0 mag in the pre-shock region (Section 4.2). With Equations (1) and (4), the value of the mass accretion rate derived from the dereddened luminosity is ${\dot{M}}_{{\rm{P}}}\gtrsim 5\times {10}^{-7}$ M_Jup yr⁻¹ (left panel in Figure 3). The PDS 70c mass accretion rate was revised to ≳1 × 10⁻⁷ M_Jup yr⁻¹ with A_Hα > 1.1 mag (middle panel in Figure 3). Note that our estimation of the mass accretion rate in this paper is limited by the detection limit of the Hβ flux, and the true mass accretion rate should be higher than the current values. In particular, we speculate that the intrinsic mass accretion rate of PDS 70c is larger than that of PDS 70b because the infrared color of K1–L' in PDS 70c with ∼2.2 mag is redder than that of PDS 70b with ∼1.1 mag (Haffert et al. 2019), resulting in a much larger value of A_Hα for PDS 70c. To better constrain the mass accretion rate of PDS 70b and c, deep multiple line observations with less extinction, such as Paβ (1.282 μm) and Brγ (2.166 μm), are preferable.

Our analysis suggests that the mass accretion rates of PDS 70b and c are ≳8× and ≳2× higher than that of the primary star (∼6 × 10⁻⁸ M_Jup yr⁻¹; Haffert et al. 2019). Recently, Thanathibodee et al. (2020) applied magnetospheric accretion models to the Hα line profile of PDS 70 and derived mass accretion rates onto the star in the range of (0.6–2.2) × 10⁻⁷ M_Jup yr⁻¹. Even with these new stellar values, the mass accretion rate of PDS 70b is still higher. Since the number of accreting planets embedded in the protoplanetary disk is currently believed to be three (LkCa 15b, PDS 70b, and PDS 70c; Kraus & Ireland 2012; Sallum et al. 2015; Wagner et al. 2018; Haffert et al. 2019), it is uncertain whether this situation (i.e., a higher planetary accretion rate) is common or rare. Numerical simulations predict that the mass accretion rate of a 1 M_Jup planet can reach up to ∼90% of the gas flow from the outer disk (Lubow & D'Angelo 2006). Furthermore, Tanigawa & Tanaka (2016) showed that the mass accretion rate of a planet in the gas gap of the disk can exceed the stellar mass accretion rate in the case of a lower planetary mass and/or a higher gas scale height (see Equation (16) in Tanigawa & Tanaka 2016). Therefore, a situation similar to the PDS 70 system can occur at a certain evolution phase in other disk systems.

Additionally we compared our results with the mass accretion rates of other young planetary-mass companions (GSC 06214–0210b, GQ Lup b, DH Tau b, and SR 12 c in Bowler et al. 2011; Zhou et al. 2014; Santamaría-Miranda et al. 2018). The mass accretion rate of GQ Lup b (${\dot{M}}_{{\rm{P}}}\,\sim 5\times {10}^{-7}$ M_Jup yr⁻¹; Zhou et al. 2014) is comparable with that of PDS 70b, while the other three objects have lower values by a few orders of magnitude (${\dot{M}}_{{\rm{P}}}\sim {10}^{-9}\mbox{--}{10}^{-8}$ M_Jup yr⁻¹; Bowler et al. 2011; Zhou et al. 2014; Santamaría-Miranda et al. 2018). These results of the higher mass accretion rates of PDS 70b and GQ Lup b could be naturally explained that these two are embedded in the circumstellar disk (MacGregor et al. 2017; Keppler et al. 2019) and can be supplied with a fresh disk material from the parent disk.

Our analysis based on the accretion shock model suggests that ${f}_{{\rm{f}}}\sim {10}^{-2}\mbox{--}{10}^{-3}$ for PDS 70b and c. If the Hα emitting regions are actually accretion shock regions, the accretion flow toward the protoplanets should be significantly converged by some mechanism. The existence of accretion shocks by converging flows has been suggested for pre-main-sequence stars such as classical T-Tauri stars, where strong stellar magnetic fields guide and collimate the accretion flows toward the magnetic poles (Koenigl 1991; Hartmann et al. 2016). The very small filling factors for PDS 70b and c may imply that magnetospheric accretion operates even in these protoplanets (for recent theoretical studies, see Batygin 2018; Thanathibodee et al. 2019). However, the existence of sufficiently strong planetary magnetic fields remains unclear, and needs to be examined observationally (e.g., radio observations with a low frequency of 10s of MHz are currently being performed; Zarka et al. 2019). If there are no magnetic fields or only a weak field in PDS 70b and c, some hydrodynamic processes should be responsible for collimating the accretion flows.

When a protoplanet does not develop a magnetosphere due to it having a weak magnetic field, and its circumplanetary disk extends to the planetary surface, the boundary layer between the protoplanet and the disk is heated due to a viscous process and it can be a source of Hα emissions. Boundary layer accretion has also been considered for pre-main-sequence stars (Bertout et al. 1988). Further theoretical and observational studies are required to identify the accretion processes. For example, high spectral resolution observations to search for the inverse P Cygni profile will help investigate the possibility of magnetospheric accretion (e.g., AA Tau; Edwards et al. 1994). Observational investigations of planetary magnetic fields at radio wavelengths will give constraints on the field strength (e.g., the Square Kilometre Array, SKA; Zarka et al. 2015). Multi-epoch observations of accreting signatures at optical/near-infrared wavelengths (e.g., VLT/MUSE or Keck/OSIRIS) will provide the information on time variability; magnetospheric accretion in pre-main-sequence stars commonly show time variability (e.g., Cody & Hillenbrand 2014).