Accretion Properties of PDS 70b with MUSE

We report a new evaluation of the accretion properties of PDS 70b obtained with the Very Large Telescope/Multi Unit Spectroscopic Explorer. The main difference from the previous studies of Haffert et al. and Aoyama & Ikoma is in the mass accretion rate. Simultaneous multiple line observations, such as Hα and Hβ, can better constrain the physical properties of an accreting planet. While we clearly detected Hα emissions from PDS 70b, no Hβ emissions were detected. We estimate the line flux of Hβ with a 3σ upper limit to be 2.3 × 10−16 erg s−1 cm−2. The flux ratio FHβ/FHα for PDS 70b is <0.28. Numerical investigations by Aoyama et al. suggest that FHβ/FHα should be close to unity if the extinction is negligible. We attribute the reduction of the flux ratio to the extinction, and estimate the extinction of Hα (AHα) for PDS 70b to be >2.0 mag using the interstellar extinction value. By combining with the Hα linewidth and the dereddening line luminosity of Hα, we derive the PDS 70b mass accretion rate to be ≳5 × 10−7 MJup yr−1. The PDS 70b mass accretion rate is an order of magnitude larger than that of PDS 70. We found that the filling factor ff (the fractional area of the planetary surface emitting Hα) is ≳0.01, which is similar to the typical stellar value. The small value of ff indicates that the Hα emitting areas are localized at the surface of PDS 70b.


Introduction
Gas giant planets growing in a protoplanetary disk gain their mass via mass accretion from the parent disk until their host star loses its gas disk (e.g., Hayashi et al. 1985). A part of the gas flow from the outer disk accretes onto gas giant planets while the rest of the flow passes over the planets toward the central star. The relationship between the host star and planets in the mass accretion rate depends on the planetary mass and the disk's properties (e.g., Lubow et al. 1999;Tanigawa & Tanaka 2016). Numerical simulations of disk-planet interactions show that the mass accretion rate of a 1 M Jup planet can reach up to ∼90% of the rate of mass accretion from the outer disk, i.e., the planetary-mass accretion rate is an order of magnitude larger than the stellar rate (Lubow & D'Angelo 2006). As predicted by planet population synthesis models (e.g., Ida & Lin 2004), when the gas accretion onto planets is sufficiently large, gas giant planets can form more than sub-Jovian planets. On the other hand, microlensing observations show that the number of sub-Jovian planets is dominant in the planetary-mass function (Suzuki et al. 2018). The mass accretion onto a planet determines the final mass of the planet (e.g., Tanigawa & Ikoma 2007) and thus investigating the planetary-mass accretion process informs our understanding of planet formation.
Direct imaging allows the study of gas accretion in young planets via H I recombination lines such as Hα at 656.28nm (e.g., Zhou et al. 2014;Sallum et al. 2015;Santamaría-Miranda et al. 2018;Wagner et al. 2018;Haffert et al. 2019). Although the gas temperature in the surface layers of protoplanetary disks is expected to be high enough (T∼10 4 K) to emit H I recombination lines (e.g., Kamp & Dullemond 2004), their emissions may be too weak to be detected due to the low density of the disk surface. Therefore, the observed Hα emissions from point-like sources reported in the literature have been attributed to the abundant cold gas around planets in the disk midplane being heated by gas accretion onto planets. In the case of TTauri stars, two important heating mechanisms in the standard accretion processes have been investigated to understand stellar accretion processes (Bertout et al. 1988;Koenigl 1991;Popham et al. 1993;Calvet & Gullbring 1998). The first is viscous heating in the boundary layer where the faster rotating Keplerian disk connects with the slowly rotating star. The other is accretion shock heating by magnetospheric accretion, in which the gas flow from the midplane of a circumstellar disk shocks the central star by a strong stellar magnetic field. Theoretical studies on planetary accretion processes have also investigated accretion shock heating at the surface of planets (e.g., Tanigawa et al. 2012;Zhu 2015;Aoyama et al. 2018).
Recently, Aoyama & Ikoma (2019) demonstrated that the Hα spectral linewidth and luminosity can be used to constrain the planetary mass and the rate of mass accretion onto protoplanets by applying radiation-hydrodynamic models of the shock-heated accretion flow (Aoyama et al. 2018). A good example of robustly detected accreting planets at tens of astronomical units from a central star is the PDS70 planetary system. This system includes a weakly accreting young star (spectral type of K7, mass of 0.8 M e , accretion rate of 6 × 10 −11 M e yr −1 , age of 5 Myr, distance of 113 pc; Pecaut & Mamajek 2016;Gaia Collaboration et al. 2018;Müller et al. 2018;Haffert et al. 2019;Keppler et al. 2019) associated with two planetary-mass companions Müller et al. 2018;Wagner et al. 2018;Christiaens et al. 2019aChristiaens et al. , 2019bHaffert et al. 2019;Isella et al. 2019). Aoyama & Ikoma (2019) estimated that PDS70b's mass and mass accretion rate are ∼12M Jup and ∼4×10 −8 M Jup yr −1 , respectively, while the values for PDS70c are ∼10M Jup and ∼1×10 −8 M Jup yr −1 . To convert the Hα luminosity into the mass accretion rate, we need to constrain the fractional area of the planetary surface emitting Hα, which is often termed the filling factor f f , and the degree of extinction A Hα of Hα. However, these values have been poorly measured in observations.
In this paper, we revisit the accretion properties of PDS70b by re-analyzing archive data used in Haffert et al. (2019). We aim to estimate the physical quantities related to planetary-mass accretion with the help of the theoretical study by Aoyama & Ikoma (2019). We will demonstrate that by combining a 3σ upper limit in the flux ratio of F Hβ /F Hα , the Hα line luminosity, and the Hα linewidth, we can constrain the filling factor f f , the Hα extinction A Hα , and the mass accretion rate onto the planet.

Archival Data
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 0 70-0 81, 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 7 42 × 7 43, a spatial sampling of 25 × 25 mas 2 , a spectral range of 4650-9300 Å, a resolving power of 2484 ± 5 at 6500.0 Å) are publicly available on the ESO archive. 13 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 spectrophotometric 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 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Å.

Stellar Spectral Subtraction
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 −20 erg s −1 cm −2 Å −1 ) 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.

Flux Measurement
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 Notes. a Flux calibration (i.e., aperture and a flux correction described in Section 2.3) is applied to values, while values are not dereddened with A Hα and A Hβ derived in Section 4.2. Dereddened Hα fluxes are >50.6 and >8.5×10 −16 ergs −1 cm −2 for PDS70b and c, respectively (see details in Section 4.2). b Line center relative to the rest wavelength of Hα at 6562.8Åis measured by Gaussian fitting. For PDS70, the line center and the 10% and 50% linewidths are −28±7, 283±28kms −1 , and 151±15kms −1 , respectively. c 3σ upper limit.
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°w ith 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. Table 1 show our re-analysis results. They are largely similar with those in Haffert et al. (2019). As already reported in Haffert et al. (2019), we also clearly detected Hα emissions from PDS70b and c, while Hβ emissions were not detected. The peak S/N of ∼20 in PDS70b (Figure 1(c)) is roughly two times better than that in Haffert et al. (2019). The Hα line fluxes with flux calibration described in Section 2.3 are 8.1±0.3×10 −16 and 3.1±0.3×10 −16 ergs −1 cm −2 for PDS70b and c, respectively. The 50% linewidths (FWHM) of ∼110kms −1 in the two planets are comparable with a spectral resolution of ∼120kms −1 in MUSE, and thus, these values could be upper limits. The main differences in our results from those in Haffert et al. (2019) are the Hα line profiles. Figure 2 shows the normalized Hα line profiles of the primary star, PDS70b, and PDS70c. All three spectra exhibit a blueshifted single Gaussian profile at the line center of −20 to −30kms −1 (−0.4 to −0.7 Å). According to the MUSE manual, this blueshift could be due to wavelength calibration errors, as the wavelength calibration errors are more than 0.4Åbased on only the arc frames taken in the morning calibrations. Note that the radial velocity of PDS70 (∼3 km s −1 : 0.07 Å; Gaia Collaboration et al. 2018) and the Keplerian velocities of PDS70b and c (∼3-4 km s −1 : 0.07-0.09 Å) are negligible.

Figure 1 and
We found a possible extended structure in the southwest of PDS70b in Figure 1(a). This structure was also reported in Haffert et al. (2019). Although the S/N for this structure is ∼3-4, the structure is exactly located at the image slicing and field splitting axis in one data cube (see Figure A1 in the Appendix). Therefore, Haffert et al. (2019) regarded this structure as an artifact. However, as described in Section 2.2, since we removed this data cube from the final image to improve the S/N of PDS70b, this structure might be marginally real. Since the study of the nature of this structure is beyond the scope of this paper, the investigations are left for future work. Furthermore, Mesa et al. (2019) recently reported the detection of a point-like feature at r∼0 12, while no significant Hα emission can be seen in our results.

Free-fall Velocity of Gas and Pre-shock Number Density of Hydrogen Nuclei
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 freefall velocity v 0 of the gas and the pre-shock number density of hydrogen nuclei n 0 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 0 and v 0 . This is because an increase of n 0 enhances the absorption of Hα emissions near the line center by shock-heated gas, while a high value of v 0 means that a flow with higher velocity passes. Thus the Hα linewidth is a function of n 0 and v 0 (see Figure 5 in Aoyama & Ikoma 2019). In Figure 3, the intersections of the Hα 10% and 50% linewidths indicate that (v 0 , n 0 ) is roughly (144 km s −1 , 3.8×10 13 cm −3 ) and (139 km s −1 , 3.2 × 10 13 cm −3 ) for PDS70b and c, respectively.

Line Flux Ratio and Extinction
The theoretical model of Aoyama & Ikoma (2019) suggests that the flux ratio of Hβ to Hα ( b a F F H H ) increases with n 0 when there is no foreground extinction, e.g., interstellar extinction (right panel in Figure 3). This is because the higher density of n 0 10 13 cm −3 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 0 or n 0 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 0 10 13 cm −3 . In contrast, a Hα absorption is negligible with the lower density of n 0 10 12 cm −3 , resulting in a smaller value of F Hβ /F Hα 0.5.
With values of n 0 and v 0 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 b F H /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.

Filling Factor, Planetary Mass, and Mass Accretion Rate
We estimate the filling factor f f using the following equation for the Hα luminosity (Aoyama & Ikoma 2019): where R p and a I H 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). a I H is a function of n 0 and v 0 , and its dependence is investigated in Aoyama et al. (2018). Using these results, we plot the lines of - 6 for PDS70b and c, respectively. To reproduce the dereddened luminosity at these cross points, the values of 10 A/2.5 /f f are estimated to be 5×10 2 and 9×10 2 , respectively. Note that the values of A Hα used here are 2.0 and 1.1mag for PDS70b and c, respectively (Section 4.2). The right panel shows the ratio of Hβ to Hα. The cross points in the left and middle panels indicate that the flux ratio is close to unity for PDS70b and c.

Planetary Mass
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 −1 . 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 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.

Origin of Planetary Extinction
The extinction of the primary star was found to be negligible ). 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 12 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 −2 , which can be translated to a molecular hydrogen surface density of ∼3×10 21 cm −2 . By using the relationship between the interstellar dust extinction and the hydrogen column density (N H /A V =1.9×10 21 cm −2 mag −1 ; 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.

Planetary Mass Accretion Rate
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 f , n 0 , and v 0 . The values of n 0 and v 0 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 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 ´- M 5 10 P 7 M Jup yr −1 (left panel in Figure 3). The PDS70c mass accretion rate was revised to 1×10 −7 M Jup yr −1 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 −8 M Jup yr −1 ; 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 −7 M Jup yr −1 . 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 ( M Ṕ . 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.

Accretion Process
Our analysis based on the accretion shock model suggests that --f 10 10 f 2 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).

Summary
We re-analyzed MUSE archive data obtained with commissioning observations, and estimated the upper limit of Hβ emissions for PDS70b. Most of the observational results are similar to those of previous studies in Haffert et al. (2019). The main difference is the planetary-mass accretion rate. We derived the accretion properties of PDS70b by assuming that the Hα emissions originate from gas accretion shock. We showed that the line flux ratio F Hβ /F Hα is useful to constrain the planetary-mass accretion rate by estimating the extinction of A Hα because the accretion rate is described as a function of A Hα . The 3σ upper limit of F Hβ /F Hα =0.28 can be translated to A Hα > 2.0mag. We then obtained a value for dereddened F Hα > 5×10 −15 ergs −1 cm −2 . With the Hα linewidth and the dereddened Hα line luminosity for PDS70b, we derived a mass accretion rate of ´- M 5 10 7 M Jup yr −1 for PDS70b. PDS70b's mass accretion rate is an order of magnitude larger than that of PDS70 with ~´-M 6 10 8 M Jup yr −1 . We also derived a filling factor f f of 0.01 for PDS70b. This result suggests that the Hα emitting areas are localized at the surface of PDS70b. Multiple line observations, especially emission lines with low extinction, such as Paβ (1.282 μm) and Paδ (2.166 μm), are useful for determining better constraints on the planetary accretion properties of young planets deeply embedded in molecular clouds or circumstellar disks. . The image size is 1″×1″. The white circles and ellipse indicate the positions of PDS70b, PDS70c, and a possible extended structure. (c) Spectra in the stripe (red) and in the off-stripe (black) at the red and black square regions in panel (a), respectively. The stripe pattern is induced by a wavelength shift possibly due to calibration errors.