Variable Accretion onto Protoplanet Host Star PDS 70

We obtained an NIR spectrum of PDS 70 using the FIRE spectrograph on the Magellan Baade Telescope at the Las Campanas Observatory in Chile. With the 06 slit, the spectrograph achieved a resolution of ∼6000 for simultaneous spectral coverage in the range 0.9–2.4 μm. We observed the star in two nodding positions (A/B), each with an exposure time of 126.8 s in the sample-up-the-ramp mode. A telluric standard star was observed immediately afterward. The average airmass and seeing of the observation were 1.1 and 10, respectively. The spectral extraction, wavelength calibration, and telluric correction were performed using the IDL-based FIRE data reduction pipeline (Simcoe et al. 2013).

PDS 70 was observed with TESS (TIC 179413040) with 2 minutes cadence for 24 days, and the photometry was reduced by the TESS pipeline. We downloaded the extracted light curve from the Mikulski Archive for Space Telescopes (MAST). The Simple Aperture Photometry (SAP) flux and the Pre-search Data Conditioning SAP flux (PDCSAP) from the pipeline appear to be the same for the star, and we choose to use the PDCSAP flux. We use the Lightkurve software (Lightkurve Collaboration et al. 2018) to extract and normalize the light curve. Figure 1 shows the 24 day light curve of PDS 70.

Zoom In Zoom Out Reset image size

We downloaded spectra of PDS 70, obtained with the High-Accuracy Radial velocity Planetary Searcher (HARPS) instrument (Mayor et al. 2003) on the European Southern Observatory (ESO) 3.6 m telescope (Program ID 098.C-0739, PI: Lagrange), from the ESO Archive Science Portal. HARPS provided 32 high-resolution (R ∼ 115,000) spectra of the star observed across ∼2 months in 2018. We used the data automatically reduced and calibrated by the HARPS pipeline for our analysis. Table 1 shows the details of the observations.

3.1.1. Rotation Period

To determine the rotation period of PDS 70, we constructed a periodogram from the TESS light curve, using the Psearch code (Saha & Vivas 2017), which combines the Lomb–Scargle and Lafler–Kinman methods. As shown in the bottom left panel of Figure 1, the periodogram shows a strong power at the period of 3.03 days. To estimate the uncertainty of the period, we fit a Gaussian to the strongest peak in the periodogram (Venuti et al. 2017). The corresponding period from the fit is 3.03 ± 0.06 days, where the uncertainty is the 1σ propagation of the fitted Gaussian width in frequency space.

We visually inspect the period by folding the light curve with the calculated period of 3.03 days. As shown in the bottom right panel of Figure 1, the period is consistent with the light curve. The amplitude of the light curve is on the order of 10%, consistent with variation due to starspots phasing in and out of the line of sight (Herbst et al. 1994; Venuti et al. 2017). Smaller dips can be seen at phase ∼0.3–0.4, suggesting that there could be other factors that modulate the light curve, such as obscuration of the inner disk warp (Bouvier et al. 2003). The 3.03 days period is also consistent with a typical rotation period of T Tauri stars (e.g., Karim et al. 2016). Therefore we conclude that the period of the light curve is the rotation period of the star since other types of variations would modulate the light curve in different timescales (e.g., Siwak et al. 2018).

We calculate the corotation radius, at which the Keplerian orbital period is equal to the stellar rotation period, and outside of which accretion cannot occur (Hartmann et al. 2016). Using the stellar mass ${M}_{\star }=0.76\pm 0.02\,{M}_{\odot }$ and radius ${R}_{\star }=1.26\,\pm 0.15\,{R}_{\odot }$ (Müller et al. 2018), we obtain the corotation radius as

$$\begin{eqnarray}&&{R}_{\mathrm{co}}={\left(\displaystyle \frac{{{GM}}_{\star }{P}^{2}}{4{\pi }^{2}}\right)}^{1/3}=6.4\pm 0.8\,{R}_{\star },\end{eqnarray} \tag{ 1 }$$

where the uncertainty is calculated with the standard error propagation.

3.1.2. Rotational Velocity

We used the Fourier method (Carroll 1933) to calculate the projected rotational velocity ($v\sin i$) of the star. The method requires isolated photospheric lines at sufficient S/N at a high spectral resolution to have a reliable line shape (Simón-Díaz & Herrero 2007). Since none of the HARPS spectra has enough S/N for the analysis, coadding the spectra is required. We first calculated the radial velocity (RV) by cross-correlating each spectrum with a PHOENIX photospheric template (Husser et al. 2013) and fitting a Gaussian to the cross-correlation function. The median and the standard deviation of the RV is $+6.0\,\pm 1.5\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$. We corrected all 32 spectra with $\mathrm{RV}=6.0\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and combined them by averaging the spectra weighted by the median S/N (see Table 1). The resulting spectrum has a median S/N of ∼46. We then selected 20 photospheric lines that are clearly isolated, including the Li i λ6708 line, and calculated a Fourier power spectrum for each line. From the first zero in the Fourier spectrum, we calculated the v sin i, adopting a limb darkening coefficient $\epsilon =0.6$, which is appropriate for spectra in the optical range (Claret 2000). Based on the median and the standard deviation of v sin i measured from these 20 lines, the rotational velocity of the star is $v\sin (i)=16.0\pm 0.5\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

We combined our measurement of v sin i with the rotation period P from Section 3.1.1 and ${R}_{\star }=1.26\pm 0.15$ R_⊙ (Müller et al. 2018), to obtain the inclination of the star as

$$\begin{eqnarray}&&i={\sin }^{-1}\left(\displaystyle \frac{{Pv}\sin (i)}{2\pi {R}_{\star }}\right)=50\pm 8^\circ \,(1\sigma ).\end{eqnarray} \tag{ 2 }$$

Keppler et al. (2019) found i = 517 ± 01 and i = 521 ± 01 for the inclination of the protoplanetary disk in two different Atacama Large Millimeter/submillimeter Array (ALMA) observations. Therefore, our results suggest that the stellar and disk rotation axes are parallel to each other within the measurement errors.

Optical hydrogen lines such as Hα have been used to determine if a T Tauri star is still accreting. However, Hα may fail at very low accretion rates due to a comparatively significant contribution from chromospheric emission (Manara et al. 2017), especially when observed in low or moderate spectral resolution. On the other hand, the He i λ10830 line has been found to be very sensitive to accretion and thus to be a good tracer for low accretors (Thanathibodee et al. 2018). In particular, the presence of a redshifted absorption component in the line profile is a direct indicator of accretion. Here we use this line to probe the state of accretion of PDS 70.

The He i λ10830 line is shown in Figure 2.⁷ In general, the He i λ10830 line is conspicuous in accreting stars. However, for low accretors observed at a medium or low spectral resolution, the Si i at 10830 Å (vacuum) could contaminate the He i line at 10833 Å. Therefore, we constructed a photospheric template of the star from interpolating in a PHOENIX model spectrum (Husser et al. 2013) with the same effective temperature and gravity as PDS 70, and convolved it first to the rotational velocity of the star and then to the resolution of the FIRE spectrograph. As shown in the left panel of Figure 2, the contribution from the Si i to the overall absorption of the He i line is negligible. The right panel of Figure 2 shows the line after subtracting the photospheric template. The line shows strong and conspicuous redshifted and blueshifted absorption features. The presence of the redshifted absorption feature is a definitive indication that the star is accreting. As with many low accretors with this type of profile, the He i emission is weak or undetectable (T. Thanathibodee et al. 2020, in preparation).

Zoom In Zoom Out Reset image size

At very low levels of accretion, the most reliable way to measure the mass accretion rate is to model the resolved profile of emission lines (Thanathibodee et al. 2019b). High resolution is needed to distinguish the chromospheric feature of the line, which appears as a narrow and mostly symmetric feature at the line center, and magnetospheric features that extend out to the star's freefall velocity. Here we use the magnetospheric accretion model of Muzerolle et al. (2001) to model the Hα line profiles from the HARPS spectra.

3.3.1. Magnetospheric Accretion Model

The physics of the magnetospheric accretion models are given in detail in Hartmann et al. (1994) and Muzerolle et al. (1998, 2001). Here we describe the basic assumptions. The model assumes that the accreting material flows along the magnetic field of the star, taken as dipolar. It assumes that the stellar rotation axis, the dipolar magnetic axis, and the disk Keplerian rotation axis are aligned. As a result, the accretion flow is axisymmetric. The free parameters of the model are the disk truncation radius (R_i), the radial width of the accretion flow on the equatorial plane (W_r), and the maximum temperature of the gas in the flow (T_max). The density in the flow is set by the mass accretion rate, $\dot{M}$, and the geometry. We solve a 16-level hydrogen atom, in which the mean intensities for the radiative rates are calculated with the extended Sobolev approximation, to obtain level populations and source functions. The line profile is calculated using a ray-by-ray method for a given inclination (i) between the magnetic axis and the line of sight.

We calculate a large grid of models varying the model parameters, as shown in Table 2. Although we have a measurement of the inclination of the rotation axis (see Section 3.1.2), we probe a larger range of inclinations to verify the axisymmetric assumption of the model; differences between the model inclination and the rotational inclination would suggest a possible misalignment between the rotation axis and the magnetic axis. We only calculate models in which the outer radius of the flow, R_i + W_r, is inside the corotation radius of the star ${R}_{\mathrm{co}}=6.4\,{R}_{\star }$. The range of $\dot{M}$ is selected to cover a typical range at which Hα starts to fail as an accretion diagnostics (few $\times {10}^{-10}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$; Thanathibodee et al. 2018). The lower end of the $\dot{M}$ range is chosen from our pregrid calculations, in which the Hα profiles can still be in emission. The range of the flow temperatures is selected to be consistent with the expected T_max at low $\dot{M}$ (Muzerolle et al. 2001; Thanathibodee et al. 2019b). In total, we calculate 64,800 models.

Table 2.  Range of Model Parameters

Parameters	Minimum	Maximum	Step
$\dot{M}$ (${10}^{-10}{M}_{\odot }\,{\mathrm{yr}}^{-1}$)	0.2	4.5	0.1, 0.5
Tmax (K)	10,000	12,000	250
Ri (R⋆)	2.0	6.0	0.4
Wr (R⋆)	0.2	2.0	0.4
i	30	75	5

Download table as:  ASCIITypeset image

3.3.2. Fitting Hα Line Profiles

The chromospheric contribution to the hydrogen lines becomes significant in low accretors. In addition, photospheric absorption lines can affect the shape of the redshifted absorption features in the line. Therefore, the photospheric and chromospheric contributions to the line need to be taken into account before modeling the line profile.

We constructed a photospheric template of the star using its normalized spectrum at the most quiescent state, during which the Hα line is symmetric and purely in emission. We replaced the Hα emission feature and any small features within 30% of the standard deviation of the flux with ${F}_{\lambda }={F}_{\lambda ,\mathrm{norm}}=1$. We then used a box filter to smooth the spectrum, resulting in a photospheric template. The left panel of Figure 3 shows a spectrum of the star and the photospheric template derived from the stellar spectra. In comparison, we plot a photospheric template interpolated from the PHOENIX model spectra (Husser et al. 2013), in the same spectral resolution and rotational velocity. The template is similar to the PHOENIX spectra, and it better reproduces the Fe i absorption line at 6569.2 Å.

Zoom In Zoom Out Reset image size

We fit each observation by calculating the rms error (RMSE) for all profiles in the grid of models. The RMSE is given by

$$\begin{eqnarray}&&\mathrm{RMSE}=\sqrt{\displaystyle \frac{\sum _{i=1}^{N}{\left({F}_{\mathrm{obs},i}-{F}_{\mathrm{model},i}\right)}^{2}}{N}},\end{eqnarray} \tag{ 3 }$$

where F_is are normalized fluxes at any given pixel and N is the total numbers of pixels in the relevant velocity range. This statistic avoids giving weight to any particular feature of the observed profile, unlike the χ² statistics, which is biased toward emission features if the deviation is normalized by the observed flux at a given pixel. We only consider the velocity range of −250 to $+400\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$, comparable to the star's freefall velocity, and exclude the region $\pm 60\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$ to avoid fitting the chromospheric feature of the line. The best-fit models are the models with the smallest RMSE, and the mass accretion rate and the accretion geometry is inferred from the weighted mean parameters of the model with RMSE ≤ 0.1. To verify that the line center is Gaussian, indicating that it arises in the chromosphere, we fit a Gaussian profile to the residual of the best fit, and then add the Gaussian profile to the photospheric profile. The right panel of Figure 3 shows an example of the best fit for one of the observed profiles. The observed line profile is reproduced by adding the best-fit magnetospheric profile with a Gaussian profile.⁸ The mass accretion rate of the star based on this observation is $1.0\times {10}^{-10}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$.

3.3.3. Accretion Variability

To explore the variability of the mass accretion rate and accretion geometry, we grouped the 32 HARPS spectroscopic observations by rotation phases, calculated from the period in Section 3.1.1, and stacked them, resulting in eight spectra representing observations in different phases. The phase ϕ is defined such that ϕ = 0 at Modified Julian Date 0.0. The grouping was done such that observations on the same night fell on one group, and the profiles were similar. As shown in Figure 4, the star shows rotational variability in the line profiles. From phase 0.00, the redshifted absorption component becomes stronger and is strongest at phase 0.49. On the other hand, the magnetospheric component of the line from phase 0.57 to 0.90 is mostly in emission, with weak to no redshifted absorption.

Zoom In Zoom Out Reset image size

To quantify the phase variation of properties of the accretion flow, we fitted the stacked spectra at each phase with the grid of models, as described in 3.3.2. The results are shown in Table 3, in which the last column shows the number of models within the RMSE ≤ 0.1 criterion. The last row of the table shows the phase-averaged values for each model parameters calculated by weighting each observation with a phase duration ${\rm{\Delta }}\phi$, where

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{i}=\displaystyle \frac{{\phi }_{i+1}-{\phi }_{i}}{2}+\displaystyle \frac{{\phi }_{i}-{\phi }_{i-1}}{2}.\end{eqnarray} \tag{ 4 }$$

Except for the flow temperature T_max, the model parameters show various degrees of variability. In particular, the observed $\dot{M}$ can vary by more than a factor of three during the rotation phase.

Table 3.  Results of the Magnetospheric Accretion Model

Phase	$\dot{M}$	Ri	Wr	Tmax	i	Number
	(${10}^{-10}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$)	(R⋆)	(R⋆)	(104 K)	(deg)	of models
0.00	1.3 ± 1.1	3.6 ± 1.0	0.5 ± 0.4	1.09 ± 0.07	47 ± 14	1229
0.23	1.7 ± 1.3	3.7 ± 1.0	0.7 ± 0.5	1.09 ± 0.07	42 ± 12	1548
0.33	1.8 ± 1.3	3.5 ± 0.9	0.7 ± 0.4	1.09 ± 0.07	50 ± 13	1441
0.49	2.2 ± 1.3	3.9 ± 0.9	0.8 ± 0.4	1.10 ± 0.07	48 ± 12	1333
0.57	0.6 ± 0.2	3.1 ± 1.3	0.2 ± 0.0	1.09 ± 0.06	48 ± 17	259
0.65	0.9 ± 0.7	2.7 ± 0.9	0.3 ± 0.2	1.10 ± 0.06	55 ± 15	432
0.84	1.0 ± 0.8	3.0 ± 0.9	0.3 ± 0.2	1.09 ± 0.06	58 ± 14	483
0.90	0.8 ± 0.7	3.0 ± 1.0	0.3 ± 0.2	1.09 ± 0.06	56 ± 15	571
Average	1.3 ± 0.5	3.4 ± 0.4	0.5 ± 0.2	1.09 ± 0.00	50 ± 6

Download table as:  ASCIITypeset image

Figure 5 shows a representative best fit for each of the observed profiles. As shown in the figure, our best-fit models show good agreement with the observed line profiles. The chromospheric fits to the residual (dashed orange lines) are centered at $v\sim 0\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$, as expected, except at ϕ = 0.57. In this case, the magnetospheric model could not well reproduce the symmetric feature of the observed profile and tends to underpredict the emission on the red wing, causing the residual (from which the chromosphere is fitted) to be redshifted. Since redshifted absorption suggests material flowing along the line of sight, the symmetric line profile suggests that there is a region of the magnetosphere with little to no flow, which is passing in front of the star at the given phase. The modifications to the current model required to test this hypothesis are beyond the scope of this paper.

Zoom In Zoom Out Reset image size

We compare the best-fit parameter values and the properties of the observed line profiles as a function of phase in Figure 6. The top two rows show the velocity at the lowest flux and the EW of the redshifted absorption feature of the Hα line, respectively. While the flow temperature is constant, general trends can be seen for other model parameters. When the absorption is weak with the smallest EW, the mass accretion rate is low, the truncation radius (R_i) is small, and the magnetosphere is thin. Thicker and larger magnetospheres correspond to stronger absorption features and higher mass accretion rates. The model inclination i, which probes the inclination of the magnetic axis, seems to vary slightly with a similar trend seen in the redshifted velocity. The slight variation in the redshifted velocity, and to some extent the inclination, seems to suggest that there could be a small misalignment between the magnetic axis and the rotation axis, resulting in a nonaxisymmetric accretion flow. Numerical simulation of accreting magnetized star with a misaligned axis (Romanova et al. 2003) suggests that even a small misalignment (${\rm{\Delta }}i\sim 5^\circ$) results in nonaxisymmetric flow, with mass preferentially flowing along a particular path. We note, however, that the uncertainty in the inclination is large compared to the model variation. Therefore, this alone, without the variability of other observed and model parameters, would not suggest the misalignment.

Zoom In Zoom Out Reset image size

The variation of the model parameters suggests that we may be probing different portions of this asymmetric flow. However, as our current model does not allow misalignment between the stellar rotation axis and the magnetic axis, we caution that such effects have to be further investigated in the future. Nevertheless, our model still provides an estimate of the mass accretion rate and its possible variation. Future spectropolarimetric observations would provide insight into the geometry of the magnetosphere (Donati et al. 2007, 2011).

We also note that in Figure 5, the chromospheric components are varying slightly between phases, especially the strength of the line, while the width and the line center are constant. We defer a discussion about chromospheric component to future studies.

Our variability analysis indicates that the mass accretion rate onto PDS 70 is in the range $(0.6-2.2)\times {10}^{-10}{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Table 3). Material accreting at this rate is slowed down through an accretion shock at the stellar surface before merging with stellar material. Emission from this shock is shown more conspicuously in UV (Calvet & Gullbring 1998; Hartmann et al. 2016).

We show optical fluxes for PDS 70 from Gregorio-Hetem & Hetem (2002) and fluxes in the Swift/Ultraviolet and Optical Telescope (UVOT) U, uvw1, and uvw2 bands (Joyce et al. 2019) in Figure 7. Except for one case, observations at uvw1 and uvw2 were obtained at different epochs. The V–I colors are similar to those of a K7 star from Pecaut & Mamajek (2013), indicating no extinction.

Zoom In Zoom Out Reset image size

We also show in Figure 7 the total flux expected for PDS 70 when the emission from accretion shocks with $\dot{M}$ between $0.8\times {10}^{-10}$ and $2.2\times {10}^{-10}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ are added to the stellar flux. The accretion shock emission is calculated using the methods of Calvet & Gullbring (1998) for the PSD 70 mass and radius, and for one characteristic value of the energy flux in the accretion column ${ \mathcal F }={10}^{12}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$ (Ingleby et al. 2013). The intrinsic stellar spectrum is taken from the WTTS HBC 427 (Ingleby et al. 2013), scaled to the stellar radius and distance. The WTTS spectrum includes emission from the stellar chromosphere, which is strong in young stars (Ingleby et al. 2011). As shown in Figure 7, the observed Swift fluxes are consistent with the accretion rates estimated from the magnetospheric modeling of Hα (Table 3).

Joyce et al. (2019) obtained an accretion rate $\dot{M}\sim 6\,\times {10}^{-12}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ using the correlation of the flux excess in the U band and the accretion luminosity from Venuti et al. (2014). As can be seen in Figure 7, the flux at the U band is dominated by stellar photospheric and chromospheric emission, as expected for low accretors (Ingleby et al. 2011), so we argue that this band is not optimal for obtaining accretion shock emission in this star. In general, relationships between the U excess and the accretion luminosity cannot be calibrated at low levels of accretion and should not be used.

Our analysis suggests that PDS 70 is accreting at a moderate rate of $\sim {10}^{-10}\,{M}_{\odot }{\mathrm{yr}}^{-1}$. Where does this mass come from? Does the outer disk gas flow across the cavity, or is there a mass reservoir to sustain the stellar accretion in the inner disk?

Numerical simulations have shown that the outer disk gas, beyond the gap opened by planets, can flow into the inner regions of the disk. The mass flow rate depends upon various factors, including the number of planets responsible for the gap, the planets' masses, the planets' accretion efficiency, and protoplanetary disk thermodynamics (Lubow & D'Angelo 2006; Zhu et al. 2011; Müller & Kley 2013). In general, these numerical studies show that the mass flow rate across the gap is 1%–100% of the mass accretion rate beyond the gap (Lubow & D'Angelo 2006; Zhu et al. 2011; Müller & Kley 2013).

In order to examine if the outer disk gas can explain PDS 70's accretion rate, we first calculate the accretion rate of an unperturbed viscous accretion disk. We adopt the disk surface density and temperature profiles used in Bae et al. (2019):

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=2.7\,{\rm{g}}\,{\mathrm{cm}}^{-2}{\left(\displaystyle \frac{R}{40\mathrm{au}}\right)}^{-1}\exp \left(-\displaystyle \frac{R}{40\,\mathrm{au}}\right),\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&T(R)=38\,{\rm{K}}{\left(\displaystyle \frac{R}{40\mathrm{au}}\right)}^{-0.24}.\end{eqnarray} \tag{ 6 }$$

In Figure 8, we present the disk accretion rate, calculated as ${\dot{M}}_{d}=3\pi \nu {\rm{\Sigma }}$ where ν is the disk kinematic viscosity defined as $\nu =\alpha {c}_{s}^{2}/{\rm{\Omega }}$. Here, α is the Shakura–Sunyaev viscosity parameter characterizing the mass transport efficiency (Shakura & Sunyaev 1973), c_s is the disk sound speed, and Ω is the angular frequency. At 70 au, beyond the common gap opened by PDS 70b and c, the accretion rate is

$$\begin{eqnarray}&&{\dot{M}}_{d}=1.5\times {10}^{-10}{M}_{\odot }{\mathrm{yr}}^{-1}\left(\displaystyle \frac{\alpha }{{10}^{-3}}\right).\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

If the disk accretion is efficient (i.e., $\alpha \gtrsim {10}^{-3}$) so the outer disk supplies gas at $\gtrsim {10}^{-10}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, it is thus possible that the stellar accretion is sustained by the mass reservoir in the outer disk, as the two planets within the gap are known to take only a small fraction of the supply ($\sim {10}^{-11}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$; Wagner et al. 2018; Aoyama & Ikoma 2019; Haffert et al. 2019; Thanathibodee et al. 2019a). However, when the mass flow rate across the gap is significantly reduced and/or the disk has a low mass transport efficiency (i.e., $\alpha \ll {10}^{-3}$), the stellar accretion rate of $\sim {10}^{-10}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ is difficult to explain with the outer disk mass reservoir. In this case, we may need to invoke an inner disk mass reservoir (e.g., a dead zone) that can feed the star for a prolonged period of time with a low efficiency (Hartmann & Bae 2018). The presence of compact submillimeter continuum emission shown in ALMA observations (Long et al. 2018; Keppler et al. 2019) may support this inner reservoir scenario.

Manara et al. (2019) used the planet population synthesis models of Mordasini et al. (2009, 2012) to make predictions for stellar mass accretion rates and disk masses and compare them with observations of the Lupus and Chamaeleon star-forming regions. Their disk models assume a viscosity parameter of $\alpha =2\times {10}^{-3}$, and planets accrete the disk gas at a fraction of the disk viscous accretion rate. With the assumed viscosity parameter, the stellar accretion could be sustained by the outer disk reservoir, subject to a decrease in the presence of planets.

Their models, however, produce a larger fraction of weak accretors than observed in transition disks. As they already pointed out, one possible explanation to this conflict is that their prescription of gas accretion onto planets overpredicts the real accretion rate, and thus reduces the mass flow rate across the gaps more than it actually would. Instead, as we suggested above, this could be reconciled with an inner disk reservoir with a small α, as the stellar accretion in this case would be less sensitive to the formation of giant planets in the outer disk. It will be interesting to run low-viscosity counterparts of planet population synthesis models to explore this possibility.

In summary, since our current understanding of protoplanetary disk accretion physics is incomplete, we cannot conclude whether PDS 70's accretion rate is sustained by an inner or an outer disk reservoir. We note that our calculations here are based on an assumption that the disk has a uniform α; other possibilities will be explored in future work. Future observations that can characterize the inner disk properties and search for potential inner disk winds, together with those that can constrain the level of turbulence in the outer disk, will help better understand the origin of PDS 70 accretion. In addition, observational searches of low accretors and accurate determination of $\dot{M}$ (Thanathibodee et al. 2018, 2019a) will help find the observational low limit of $\dot{M}$ to compare with expectations of exoplanets' population models.

The presence of a subcontinuum blueshifted absorption feature in the He i λ10830 line has been attributed to winds (Edwards et al. 2003, 2006). In particular, the narrow blueshifted absorption is interpreted as a wind coming from the inner disk. The coexistence of blueshifted and redshifted absorption with very weak emission in the line center has been observed for a few stars in the high-resolution survey of the He i λ10830 line (Edwards et al. 2006). However, most of the stars in the surveys, and all of the stars with blue+red absorption, have high levels of accretion, $\dot{M}\geqslant {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Interestingly, PDS 70, with its very low accretion level, also shows a similar type of profile. In fact, in our ongoing survey of very low accretors, more than ∼10% of disk-bearing T Tauri stars with weak Hα show blue+redshifted absorption. The velocity of the blueshifted absorption ($\sim -85\,\,\mathrm{km}\,{{\rm{s}}}^{-1}$) is consistent with the wind coming from the inner disk.

We can estimate the mass-loss rate from the wind by adopting the procedure outlined by Calvet (1997). For a line formed in the wind, the optical depth at a given velocity v is given by

$$\begin{eqnarray}&&\tau =\displaystyle \frac{\pi {e}^{2}}{{m}_{e}c}\displaystyle \frac{{fc}}{{\nu }_{0}}\displaystyle \frac{{n}_{l}(v)}{{dv}/{dz}},\end{eqnarray} \tag{ 8 }$$

where f is the oscillator strength, ν₀ is the line frequency, n_l is the number density of the lower level, and dv/dz is the velocity gradient. The mass-loss rate is given by

$$\begin{eqnarray}&&{\dot{M}}_{w}\sim {\rm{\Delta }}{Av}\mu {m}_{{\rm{H}}}{n}_{{\rm{H}}}(v)\sim {\rm{\Delta }}{Av}\mu {m}_{{\rm{H}}}\eta {n}_{l}(v),\end{eqnarray} \tag{ 9 }$$

where ΔA is the cross-sectional area of the wind at v, μ is the mean molecular weight, n_H is the number density of hydrogen, and $\eta \equiv {n}_{{\rm{H}}}/{n}_{l}$. To the first approximation, ${dv}/{dz}\sim v/{R}_{\star }$ and ${\rm{\Delta }}A\sim \pi {(2{R}_{\star })}^{2}$ (Calvet 1997). We estimate the parameter η by calculating a typical fraction between the number density of the lower level of He i λ10830 to the total hydrogen number density. We use the C17.01 release of the software Cloudy (Ferland et al. 2017) to calculate the level populations in a slab of gas assuming that the ionization radiation is an X-ray with a blackbody temperature of 5 × 10⁶ K and ${L}_{X}=5\times {10}^{29}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, consistent with the luminosity measured by Swift (Joyce et al. 2019).⁹ Assuming that the wind is ∼10 R_⋆ from the star with a thickness of 1 R_⋆, we find that η ∼ 10⁷ for a relevant range of n_H across the slab. With τ ∼ 0.5, estimated from the depth of the feature, and μ = 2.4, we find that ${n}_{l}(v=85\,\,\mathrm{km}\,{{\rm{s}}}^{-1})\sim 55\,{\mathrm{cm}}^{-3}$ and ${\dot{M}}_{w}\sim {10}^{-11}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}.$ This mass-loss rate is consistent with that expected from an MHD inner disk wind in which ${\dot{M}}_{w}\sim 0.1{\dot{M}}_{\mathrm{acc}}$, suggesting that the blueshifted absorption in He i λ10830 forms by a similar mechanism as in high accretors (Calvet 1997; Edwards et al. 2006; Kwan et al. 2007). Since the blueshifted velocity is high, it is unlikely that the feature is formed in photoevaporative winds (see Alexander et al. 2014). The features in the He i λ10830 lines of low accretors will be described in detail in a future study (T. Thanathibodee et al. 2020, in preparation).

Our detection of wind and accretion signatures confirms the existence of gas in the inner disk around PDS 70. Long et al. (2018) searched for first overtone CO lines in low-resolution near-IR spectra and could not find them, from which they inferred that the inner disk of PDS 70 was gas poor. However, the lack of detection could be due to the difficulty of separating the disk emission from the intrinsic photospheric CO absorption lines (Calvet et al. 1991). Observations of ground state CO lines or of fluorescent H₂ lines can help confirm the presence of gas in the inner disk of PDS 70.