# The Planetary Accretion Shock. III. Smoothing-free 2.5D Simulations and Calculation of Hα Emission

Gabriel-Dominique Marleau, Rolf Kuiper, William Béthune, Christoph Mordasini

► To cite this version:

Gabriel-Dominique Marleau, Rolf Kuiper, William Béthune, Christoph Mordasini. The Planetary Accretion Shock. III. Smoothing-free 2.5D Simulations and Calculation of Hα Emission. The Astrophysical Journal, 2023, 952 (1), pp.89. 10.3847/1538-4357/accf12. hal-04210292

HAL Id: hal-04210292

https://hal.science/hal-04210292v1

Submitted on 18 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# The Planetary Accretion Shock. III. Smoothing-free 2.5D Simulations and Calculation of Hα Emission

Gabriel-Dominique Marleau<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />✉, Rolf Kuiper<fnref n="2" />✉, William Béthune<fnref n="5" />, and Christoph Mordasini<fnref n="3" />✉

Received 2022 September 1; revised 2023 April 5; accepted 2023 April 19; published 2023 July 19

## Abstract

Surveys have looked for Hα emission from accreting gas giants but found very few objects. Analyses of the detections and nondetections have assumed that the entire gas flow feeding the planet is in radial freefall. However, hydrodynamical simulations suggest that this is far from reality. We calculate the Hα emission from multidimensional accretion onto a gas giant, following the gas flow from Hill sphere scales down to the circumplanetary disk (CPD) and the planetary surface. We perform azimuthally symmetric radiation hydrodynamics simulations around the planet and use modern tabulated gas and dust opacities. Crucially, contrasting with most previous simulations, we do not smooth the gravitational potential but do follow the flow down to the planetary surface, where grid cells are 0.01 Jupiter radii small. We find that roughly only 1% of the net gas inflow into the Hill sphere directly reaches the planet. As expected for ballistic infall trajectories, most of the gas falls at too large a distance on the CPD to generate Hα. Including radiation transport removes the high-velocity subsurface flow previously seen in hydrodynamics-only simulations, so that only the free planet surface and the inner regions of the CPD emit substantial Hα. Unless magnetospheric accretion, which we neglect here, additionally produces Hα, the corresponding Hα production efficiency is much smaller than usually assumed, which needs to be taken into account when analyzing (non)detection statistics.

*Unified Astronomy Thesaurus concepts: Accretion (14); Radiative transfer (1335); Hydrodynamical simulations (767); Shocks (2086); Planet formation (1241); H alpha photometry (691)*

## 1. Introduction

Over the last roughly 3% of a millennium, hundreds of extrasolar super-Jupiters have been discovered (Zhu & Dong 2021). However, only a very few of those objects are young (≤100 Myr). These are accessible mostly at large separations from their host star through direct imaging (~10–100 au; e.g., Wagner et al. 2019; Vigan et al. 2021). A few more are predicted to be on the verge of being discovered (Asensio-Torres et al. 2021), but this would not change the overall abundance qualitatively. Despite dedicated surveys (Cugno et al. 2019; Xie et al. 2020; Zurlo et al. 2020; Huélamo et al. 2022; Follette et al. 2023), even fewer gas giants have been caught accreting. It is a robust theoretical prediction that gas undergoing a shock at a velocity $v_0$ above a critical value, $v_0 > v_{H\alpha, \text{crit}} \approx 30 \text{ km s}^{-1}$, will emit hydrogen lines as the hydrogen ionized by the shock recombines and cools in the geometrically thin postshock region (Aoyama et al. 2018, 2020). There are only a few planetary-mass companions with observed line emission clearly linked to a shock: PDS 70 b and c (Wagner et al. 2018; Haffert et al. 2019) at Hα and Delorme 1(AB)b at several lines (Eriksson et al. 2020; Betti et al. 2022a, 2022b; Ringqvist et al. 2023). However, only the PDS 70 planets are found in a gas disk, and the others are effectively isolated. At newly-discovered AB Aur b, pointlike Hα emission is observed but scattering of stellar photons

cannot be excluded as the source (Currie et al. 2022; Zhou et al. 2022).

Different factors can explain the scarcity of Hα-emitting accreting planets detected at large separations. Most planets are possibly forming closer in to their star, as classically expected from core accretion (e.g., Thommes et al. 2008; see also review in Emsenhuber et al. 2021), and thus inside the inner working angle of current-generation detectors (Close 2020). Planets could also be accreting episodically, with only brief and therefore unlikely to be caught periods of detectably high accretion (e.g., Brittain et al. 2020). Another possibility is that the Hα and the other hydrogen lines are absorbed by the protoplanetary disk (PPD). Since massive planets open a deep gap in the gas and dust distributions, this might not be an important effect. As an example, when scaling the simulations of Sanchis et al. (2020) to the PDS 70 system, extinction by the PPD is negligible for their specific accretion rate (see discussion in Section 7.6 of Marleau et al. 2022). Alternatively, the gas and dust flowing onto the planet could absorb the Hα, but Marleau et al. (2022) estimated that this not important for most planet accretion rates and masses.

Two more factors are of particular importance yet underappreciated. One is that line emission in the planetary-mass regime, especially in the absence of magnetospheric accretion, is intrinsically weaker than for stars for a given mass infall rate (Aoyama et al. 2021).

The second factor is that most of the infalling and supersonic gas likely does not reach the planet directly but rather shocks onto a circumplanetary disk (CPD) at a significant fraction of the Hill radius $R_{\text{Hill}}$ away from the planet. This consequence of

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germany</footnote>

<footnote n="2">² Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1, D-47057 Duisburg, Germany</footnote>

<footnote n="3">³ Physikalisches Institut, Universität Bern, Gesellschaftsstr. 6, CH-3012 Bern, Switzerland</footnote>

<footnote n="4">⁴ Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="5">⁵ DAAA, ONERA, Université Paris Saclay, F-92322 Châtillon, France</footnote>
<table><caption>**Table 1**
Comparison of Some (Radiation) Hydrodynamical Studies of Accreting Gas Giants</caption><thead><tr><th>Hα<sup>a</sup></th><th>Study and Dimensionality</th><th>Domain Size<sup>b</sup></th><th>Thermodynamics<sup>c</sup></th><th>Smoothing<sup>d</sup> ε<sub>grav</sub></th><th>Resolution<sup>d</sup> Δr</th></tr></thead><tbody><tr><td></td><td>Machida & Kokubo (2008)</td><td>3D (30, 120, 8) H<sub>P</sub></td><td>Isothermal</td><td>0.0060 H<sub>P</sub></td><td>0.0070 H<sub>P</sub></td></tr><tr><td></td><td>Tanigawa et al. (2012)</td><td>3D (24, 24, 6) H<sub>P</sub></td><td>Isothermal</td><td>0.0007 H<sub>P</sub></td><td>0.0004 H<sub>P</sub></td></tr><tr><td></td><td>Béthune & Rafikov (2019b)</td><td>3D 128 R<sub>P</sub></td><td>Isothermal</td><td>None</td><td>0.1 R<sub>J</sub></td></tr><tr><td></td><td>Fung et al. (2019)</td><td>3D Global</td><td>Isothermal<sup>e</sup></td><td>0.05 H<sub>P</sub></td><td>0.02 R<sub>Hill</sub></td></tr><tr><td>*</td><td>Szulágyi & Ercolano (2020)</td><td>3D Global</td><td>FLD</td><td>20 R<sub>J</sub></td><td>0.9 R<sub>J</sub></td></tr><tr><td></td><td>Dong et al. (2021)</td><td>2.5D 10 R<sub>P</sub></td><td>Isothermal</td><td>None</td><td>0.001 R<sub>P</sub></td></tr><tr><td>*</td><td>Takasao et al. (2021)</td><td>2.5D 0.03 R<sub>Hill</sub> ≈ 100 R<sub>P</sub></td><td>γ = 1.01, 1.05<sup>f</sup></td><td>None</td><td>0.005 R<sub>P</sub></td></tr><tr><td>*</td><td>Maeda et al. (2022)</td><td>3D (24, 24, 6) H<sub>P</sub></td><td>Isothermal</td><td>0.0002 H<sub>P</sub></td><td>0.0004 H<sub>P</sub></td></tr><tr><td>*</td><td>This work</td><td>2.5D 1 R<sub>Hill</sub> ≈ 4100 R<sub>J</sub></td><td>FLD</td><td>None</td><td>0.001 R<sub>J</sub></td></tr></tbody></table>

**Notes.** Particularly commendable settings are highlighted in bold. All simulations reported here are for gas giants and/or have $q_{th} \gtrsim 2$ (see Section 3.3 for this parameter and others mentioned here). Dong et al. (2021) is included for comparison, even though they did not allow infall from the outer edge. We take the case $q_{th} = 4$ of Béthune & Rafikov (2019b) and Fung et al. (2019), $M_{P} = 3 M_{J}$ of Szulágyi & Ercolano (2020), and $R_{Hill}/H_{P} = 1.36$ ($q_{th} = 7.5$) of Maeda et al. (2022) because they are closest to our fiducial values (Table 2).

<sup>a</sup> An asterisk marks studies predicting Hα emission (through very different approaches).

<sup>b</sup> Three values: size in (x, y, z) in Cartesian coordinates; one value: radial extent in polar/spherical coordinates. “Global” refers to work simulating at least a significant radial ring of the PPD.

<sup>c</sup> Simulations with radiation transfer use FLD and have $\gamma = 1.43$.

<sup>d</sup> Smallest smoothing length of the gravitational potential and smallest radial cell size if they are not constant or, in nested-grid simulations, if they depend on the grid level.

<sup>e</sup> They also perform adiabatic simulations for comparison, but this is likely far from reality.

<sup>f</sup> “Nearly isothermal” would not be an appropriate term because they obtain an extremely hot postshock region (Section 5.2).

angular momentum conservation holds for matter inflowing due only to the action of gravity (i.e., ballistically; Mendoza et al. 2009). In the planet formation context, it was pointed out by Tanigawa et al. (2012), who found in their isothermal hydrodynamical simulations that the gas hitting the CPD around a Jupiter-mass planet was spread uniformly over a region of size ~0.1 $R_{Hill} \sim 100 R_{J}$. This shock should heat up the CPD and possibly make it detectable in the near-infrared (Szulágyi et al. 2019). However, for typical planet masses, the velocity of the gas at the shock $v_{0} \sim 1/\sqrt{r}$ (see Equation (2) below) is too low at large distances for significant line emission from the shock (i.e., it has $v_{0} < v_{Hα, crit}$). This implies that only the small fraction of the mass inflow that hits the planetary surface and the innermost region of the CPD would be responsible for line emission. However, these regions were not resolved in their simulations. Similarly, other recent work addressing Hα generation is limited in different ways, as we review in Section 2.

In this work, we study mass infall onto the planetary surface and the CPD, including its innermost regions. To enable an appropriately high resolution with $\Delta r \ll R_{J}$, we make the compromise of a simplified dimensionality but otherwise use methods matching or improving previous studies. We then predict the Hα emission using detailed shock emission models designed for the planetary regime.

This paper is structured as follows. Section 2 summarizes related studies in the literature, and Section 3 presents our physical model and numerical methods. In Section 4, we discuss the properties of the flow, estimate the amount of generated Hα, and assess the effect of varying some input parameters. In Section 5, we compare our results to previous work before presenting a more general discussion. Finally, in Section 6, we summarize our findings and conclude.

## 2. Previous Studies of Accreting Planets

Previous studies have looked at the accretion flow toward a forming planet. Pioneering work was presented by Ayliffe &

Bate (2009a, 2009b, 2012). However, most investigations (e.g., those and Cimerman et al. 2017; Kurokawa & Tanigawa 2018; Béthune & Rafikov 2019b; Schulik et al. 2019, 2020; Mai et al. 2020; Bailey et al. 2021; Moldenhauer et al. 2021; Krapp et al. 2022; Moldenhauer et al. 2022) were concerned with low planet masses (super-Earths to at most $1 M_{J}$), for which no Hα from an accretion shock can be expected because the infall velocity is too low (see Equation (2); Aoyama et al. 2018). The few works exploring accretion onto gas giants have limited spatial resolution ($\Delta r \sim R_{J}$) and, most problematically, a sizable gravitational potential smoothing length of $\epsilon_{grav} \sim 10 R_{J}$ or larger, which even modern computational resources impose (e.g., Machida & Kokubo 2008; Tanigawa et al. 2012; Szulágyi 2017; Fung et al. 2019; Lambrechts et al. 2019; Szulágyi & Ercolano 2020). Table 1 provides an overview. This significant $\epsilon_{grav}$ already affects the flow on scales of a few times $\epsilon_{grav}$ (i.e., out to roughly $r \sim 30-100 R_{J}$ or more) by weakening the effective mass of the planet. This could thus qualitatively change the flow pattern at the length scales that set the shock velocity, which is a sensitive factor for the Hα emission. In Table 1, studies predicting Hα emission are highlighted by an asterisk.

A notable exception to the issue of a large smoothing of the gravitational potential is the work of Takasao et al. (2021), who considered the full potential ($\epsilon_{grav} = 0$) of their $12 M_{J}$ planet. Their cell size, $\Delta r \sim 0.01 R_{J}$, close to the planet is also adequate to resolve the flow details. However, they did not include radiation transfer and adopted a heat capacity ratio of $\gamma = 1.01$ or 1.05. We will compare with their work in Section 5.2 and find crucial differences in the postshock structures and thus the emission of Hα. Therefore, including radiation transfer while setting $\epsilon_{grav} = 0$ seems desirable for more realistic simulations.

In this work, we complement the studies in the literature by considering the full nonsmoothed potential with a high spatial resolution while including radiation transport. As a compromise, instead of resolving the 3D structure of the flow, we assume axisymmetry around the planet. We describe our model in detail in the next section.
# 3. Physical Model and Numerical Methods

## 3.1. *Approach: Local Simulations in 2.5D*

We consider a super-Jupiter forming by runaway accretion in a PPD and study the flow around the planet from Hill sphere scales down to length scales much smaller than a Jupiter radius. The dynamical timescales near the planet are much shorter than the dynamical or even viscous timescales of the PPD. Therefore, we do not evolve the background PPD and take it as a fixed boundary condition, assuming a circular orbit for the planet. We expect the flow around the planet to reach a quasi-steady state over a few freefall times from the Hill sphere down to the planetary surface. While this is not a true steady state because of accretion, this state may represent an instantaneous snapshot in the accretion history of the planet. We also assume that a quasi-steady state is reached much faster than the mass and accretion rate of the planet changes while it forms.

Initially, the domain contains a negligible amount of mass and thermal energy (see Section 3.2.4). However, there is no well-defined final amount of mass and thermal energy because the quasi-steady state is an essentially constant accretion flow ($d\dot{M}/dt \sim 0$), not the absence of accretion ($\dot{M} \sim 0$). Consequently, we do not attempt to predict quantities such as the CPD thickness, size, or temperature or the interior luminosity of the accreting planet. They will likely depend on the simulation history, which is not guaranteed to be equivalent to a global calculation of planet formation in an evolving PPD. Instead, quantities such as the CPD thickness or planet luminosity should be seen as independent parameters that can be measured in the simulation results. They change slowly in the quasi-steady state.

The physical domain of our simulations extends to the Hill sphere. This is a compromise between, on the one hand, simulating the whole PPD structure, which would be computationally expensive but also sensitive to several poorly constrained modeling choices (e.g., strength and spatial dependence of viscosity, presence of disk winds, dust grain size evolution, and feedback on the gas), and, on the other hand, considering a region only several planetary radii large, which would let the chosen boundary conditions at that location determine the gas flow too strongly. Our assumption of axisymmetry around the planet will break down significantly further away than the Hill sphere, where the star's gravity dominates. Therefore, simulating out to $R_{\text{Hill}}$ is a natural choice.

This paper extends to a more realistic geometry our previous work in this series (Marleau et al. 2017, 2019; hereafter Papers I and II), in which we simulated and analyzed the properties of the accretion flow and shock with highly resolved 1D models for purely radial infall. The emphasis in the present work is on the flow geometry when including angular momentum conservation.

## 3.2. *Numerical Methods*

To solve the radiation hydrodynamics, we use the open-source (magneto)hydrodynamics code *Pluto* (Mignone et al. 2007, 2012) in combination with the radiation transport package *Makemake* (Kuiper et al. 2020). The Courant-Friedrichs-Lewy number is kept at CFL = 0.4, except for the simulation of Section 4.4, which uses CFL = 0.33. We use the nonequilibrium (two-temperature) flux-limited diffusion (FLD) module (Kuiper et al. 2010, 2020) as in Papers I and II. As argued in Tanigawa et al. (2012), we do not include viscosity,

since we study the supersonic flow, not the CPD structure. Also, self-gravity is negligible. Makemake has been extensively tested and used in a variety of contexts (see Kuiper et al. 2020 and references therein).

We do not set floor values on quantities such as the density or pressure. The only exception is a minimum on the radiation temperature $T_{\text{rad}} = (E_{\text{rad}}/a_r)^{1/4} \sim 10^{-10}$ K (where $a_r$ is the radiation constant) during the numerical iterations to solve the FLD equation system to prevent $E_{\text{rad}} \le 0$ K from ever being reached. The smallest temperatures in the converged profiles are, however, much higher ($T \gtrsim 10$ K), as expected.

### 3.2.1. *Coordinate System, Domain Size, Forces, and Resolution*

We consider a simulation region in the $r-\theta$ (poloidal, or vertical) plane centered on the planet and averaged over $\phi$, where $r, \theta$, and $\phi$ are the usual radial, spherical polar (colatitudinal), and azimuthal coordinates. We simulate the upper hemisphere and assume symmetry at the midplane ($\theta = 90^\circ$). The radial grid ranges from $r_{\text{min}} = 1.9 R_J$ to $r_{\text{max}}$, with $r_{\text{max}}$ a near-unity factor of $R_{\text{Hill}}$. Our fiducial simulation has $r_{\text{max}} = R_{\text{Hill}}$. With $r_{\text{min}} = 1.9 R_J$, the shock that defines the planet surface is usually at $R_p \approx 2 R_J$, changing only slowly during the simulation. This commonly used size is thought to be appropriate for forming or young planets (Marley et al. 2007; Mordasini et al. 2012a; Zhu 2015).

We include the gravity of the planet by adding a radial acceleration $g = -GM_p/r^2$ everywhere, where $M_p$ is the constant mass of the planet. Thus, we avoid any smoothing and set $\epsilon_{\text{grav}} = 0$. Since the mass in the whole simulation domain is always negligible, no self-gravity is needed. We do not directly include the vertical component of the star's gravity in the simulation. Above the planet, it would dominate over the planet's gravity for $z \gtrsim 3^{(1/3)} R_{\text{Hill}}$, making the star's gravity only a small correction to the dynamics in the region we simulate. However, we do impose the appropriate density stratification at the outer edge of the domain, which is sufficient (see Equation (1) below).

We follow the azimuthal component of the velocity $v_\phi$ despite the azimuthal symmetry, making the simulation 2.5D. The simulations include the Coriolis and centrifugal terms due to the planet's Keplerian orbit around the star. This is done in a linearized, conservative form (Kley 1998) by enabling the `ROTATING_FRAME` option of *Pluto*.

We use a very fine radial gridding close to the inner edge with $\Delta r = 10^{-3} R_J$, where the atmosphere is in rotation-modified hydrostatic equilibrium. The cell size increases smoothly outward. Near the shock that terminates the atmosphere and defines the radius of the planet, cells have $\Delta r \sim 10^{-2} R_J$. Beyond $r \approx 2.5 R_J$, the cell size increases logarithmically with 76 cells per decade, reaching $\Delta r \approx 100 R_J$ at $R_{\text{Hill}}$. Appendix A gives further details.

The standard polar grid is uniform, with $N_\theta = 181$ cells from pole to equator ($\Delta\theta \approx 0.5^\circ$). Simulations with $N_\theta = 91$ and even $N_\theta = 51$ ($\Delta\theta \approx 1^\circ.8$) yielded the same results overall. The only difference is that only in the middle- and high-resolution simulations is a thin supersonic flow beneath the surface of the CPD visible. We discuss this in Appendix B.

### 3.2.2. *Boundary Conditions*

The radial boundary conditions are described in the next two sections. In the polar ($\theta$) direction, we use reflective and
equatorially symmetric boundary conditions at the pole and midplane, respectively.

*Boundary conditions at the outer edge.* We assume the surrounding PPD to be vertically isothermal and in hydrostatic equilibrium. We do not include the stellar potential, but we fix the density at the outer edge to

$$\rho(r_{\text{max}}, \theta) = \rho_{\text{mid}} \exp \left[ -0.5 \left( \frac{z}{H_P} \right)^2 \right], \quad (1)$$

where $z = r \cos \theta$ is the height above the midplane, and $\rho_{\text{mid}} = \Sigma/(\sqrt{2\pi H_P})$ is the midplane volume density in the gap with constant surface density $\Sigma$ (see text below Expression (7)). Equation (1) thus correctly mimics the influence of the central star.

The poloidal components of the velocity are set as follows at $r_{\text{max}}$. We set $v_{\theta} = 0$. The radial velocity has $dv_r/dr = 0$ to allow both inflow ($v_r < 0$) and outflow ($v_r > 0$), but it is limited in magnitude on the negative side to the freefall velocity from infinity,

$$v_{\text{ff},\infty}(r) = \sqrt{\frac{2GM_{\text{p}}}{r}} = 60 \text{ km s}^{-1} \sqrt{\frac{M_{\text{p}}/2 M_{\text{J}}}{r/2 R_{\text{J}}}}. \quad (2)$$

Locally in the rotating frame, the flow of the PPD reduces to a simple linear shear (Hill 1878; Goldreich & Lynden-Bell 1965), which we average along $\phi$ to obtain the azimuthal component of the velocity. Taking $x$ to point away from the star along the star–planet direction and $y$ along the orbit of the planet, the shear is $v_y = -q\Omega_0 x$, where the Keplerian orbital angular frequency of the planet is $\Omega_0 = \sqrt{GM_*/a^3}$ for a semimajor axis $a$, and $q = 3/2$ for a Keplerian potential. Since $x = R \sin \phi$, with $R$ the cylindrical radius, we have $v_\phi = -q\Omega_0 R \sin^2 \phi$, with a $\phi$ average of $\langle v_\phi \rangle = -\frac{1}{2}q\Omega_0 R$. We therefore set

$$v_{\phi}(r_{\text{max}}) = \langle v_{\phi} \rangle = -\frac{3}{4} \Omega_0 r_{\text{max}} \sin \theta. \quad (3)$$

This is an approximation, since a shear flow is an exact description of the gas motion only close to the planet while neglecting its presence at the same time.

Averaging the radial component the same way gives $\langle v_r \rangle = 0$, which we obviously do not use for $v_r$ because it would prevent accretion into the domain. In reality, due to the planet’s gravity, there is no pure shear flow but rather complex horseshoe orbits with “U-turns” and other features that can be captured only in 3D (e.g., Tanigawa et al. 2012; Schulik et al. 2020). Thus, setting $dv_r/dr = 0$ as detailed above is a simple attempt to circumvent the limitation of a formally averaged 2.5D approach.

Finally, we take a zero-gradient boundary condition for the gas pressure ($dP/dr = 0$), which, however, is unimportant because the gas is supersonic for $r_{\text{max}} < R_{\text{Bondi}}$, which will hold for our cases of interest. We set $d(r^2 E_{\text{rad}})/dr = 0$ for the radiation energy density $E_{\text{rad}}$; if the radiation is free-streaming at $r_{\text{max}}$ (as it does turn out to be), this corresponds to a zero-gradient condition on the luminosity (Paper I).

*Boundary conditions at the inner edge.* Young planets have been observed to spin at 5%–20% of their breakup frequency (e.g., Bryan et al. 2018, 2020). Therefore, we let the planet rotate at $r_{\text{min}}$ as a solid body by setting $v_\phi(\theta) = f_\omega \omega_{\text{crit}} r_{\text{min}} \sin \theta$, where the critical or breakup frequency is given by (see, e.g.,

Section 4 of Paxton et al. 2019)

$$\omega_{\text{crit}} \approx \sqrt{\frac{GM_{\text{p}}}{R_{\text{p}}^3}}, \quad (4)$$

with the normalized spin $f_\omega$ set to 0.1. Planets spinning at near-breakup rates ($f_\omega \approx 0.8$) might shed mass more than accrete (Dong et al. 2021; Fu et al. 2023), but smaller values should not significantly influence the transport of either mass or momentum in the CPD, and certainly not in the supersonic part of the flow. Therefore, we do not vary $f_\omega$.

As in Papers I and II, the inner edge is closed, without a flow of matter. Therefore, we set $d\rho/dr = 0$ at $r_{\text{min}}$ and use a reflecting condition on the radial velocity, $v_r(r_{\text{min}}^-) = -v_r(r_{\text{min}}^+)$, where $r_{\text{min}}^{+(-)}$ is above (below) $r_{\text{min}}$. This lets an atmosphere in equilibrium build up beneath the settling zone. We use a no-slip boundary condition, $v_\theta = 0$. The pressure is determined by hydrostatic equilibrium in the presence of rotation,

$$\frac{dP}{dr} = -\rho \left( g - \frac{v_{\phi}^2}{r} \right), \quad (5)$$

where $g = -GM_{\text{p}}/r^2$ (see above).

At the inner edge, we choose a small luminosity $L(r_{\text{min}}) = 10^{-7} L_\odot$ and impose accordingly across the interface $dE_{\text{rad}}/dr = L/(4\pi r^2)/D_{\text{F}}$, where $D_{\text{F}} = \lambda c/(\kappa_{\text{R}} \rho)$ is the local diffusion coefficient, with $\lambda$ the flux limiter (see details in Paper I), and $\kappa_{\text{R}}$ is the Rosseland mean opacity. The luminosity increases outward due to the compression of the accreting gas (Paper II).

### 3.2.3. Microphysics

As in our previous work, the Rosseland and Planck mean opacities are taken from Malygin et al. (2014) for the gas and Semenov et al. (2003; model nrm_h_s) for the dust, with dust sublimation as in Isella & Natta (2005). The maximal dust-to-gas mass ratio is $f_{d/g} = 10^{-4}$. This reduction with respect to the ISM value reflects the fact that the “pressure bump” induced by the planet in the PPD could keep opacity-carrying grains out of the gap (Drażkowska et al. 2019; Chachan et al. 2021; Karlin et al. 2023; but see Szulágyi et al. 2022). This is, however, uncertain, and $f_{d/g}$ could be varied in future work.

For the equation of state (EOS), we use a perfect gas with a constant mean molecular weight $\mu = 2.3$ and adiabatic index $\gamma = 1.4$. appropriate for a solar mixture of H$_2$ and He with a hydrogen mass fraction $X = 0.75$. In Papers I and II, we showed that the choice of $\mu$ and $\gamma$ does not affect the hydrodynamic structure of the accretion flow (as seen also in Chen & Bai 2022). The CPD properties, particularly the thickness, might be especially affected by $\mu$, but the accretion history is likely as important a factor. We recall that we do not wish to predict the quantitative properties of the CPD here. Thus, the choice of the EOS will not bear qualitatively on the results.

### 3.2.4. Reaching a Quasi-steady State

Gas is free to flow into the simulation domain from the outer edge. For a steady state to be reached, at least a few global freefall times need to elapse. Evaluating the freefall time from any radius $r$ to a much smaller position at $r = r_{\text{max}}$ (e.g.,
Mungan 2009) yields the global freefall time

$$t_{\text{ff,glob}} = \pi \sqrt{\frac{r_{\text{max}}^3}{8GM_{\text{p}}}}. \qquad (6)$$

By the definition of freefall, Equation (6) ignores angular momentum. In reality, the latter will reduce the radial velocity of the gas and thus increase its fall time. For reference, if $r_{\text{max}}$ is always set to $R_{\text{Hill}}$, the freefall time becomes $t_{\text{ff,glob}} = (2\pi/\Omega_0)/\sqrt{96}$, or one-tenth of an orbital period.

We initialize the simulation with small density and temperature values that decrease outward independently of angle and set $T_{\text{rad}} = T$ (where T is the gas temperature), $v_r = v_\theta = 0$, and $v_\phi = -0.75\Omega_0 R$ throughout. This state is quickly “forgotten” over a timescale comparable to $t_{\text{ff,glob}}$ as the gas begins to fall due to the gravity of the planet. The accreting gas accumulates in a CPD whose outer edge is defined by a radial shock and grows slowly over time.

To speed up the computation, we use two phases. In Phase I, while we let the large-scale and supersonic flow reach a quasi-steady state that erases the initial conditions, we do not compute the hydrodynamics of the innermost region, whose early-time properties will be unimportant once the gas from $r_{\text{max}}$ reaches it. We use Pluto’s FLAG_INTERNAL-ARY to make the cells between $r_{\text{min}}$ and a “freeze radius” $r_{\text{frz}}$ inactive. We set $r_{\text{frz}} = 10 R_{\odot}$. Importantly, the inactive cells are ignored when determining the hydrodynamics time step<fnref n="6" /> $\Delta t$. From the Courant condition, $\Delta t$ increases with cell size and decreases with temperature. Therefore, not having to take the innermost cells into account, which are the smallest (Figure 8) and the hottest (Figure 10), speeds up the computations considerably. The radiation transport is always solved over the full domain, both in Phase I and the subsequent Phase II (described next).

After a few global freefall times in Phase I, we restart the simulation but now evolve the density and velocity everywhere as usual. This is Phase II. After a brief transition period, no features remain at $r_{\text{frz}}$, and all quantities ($\rho$, $v$, $P$, and $T_{\text{rad}}$) are smooth. We run Phase II for thousands of freefall times from $r_{\text{frz}}$ (numerical details are given in Section 3.4). By Equation (6), the freefall time at $r_{\text{frz}}$ down to $r=0$ is roughly $t_{\text{ff,frz}} = (r_{\text{frz}}/R_{\text{Hill}})^{1.5}t_{\text{ff,glob}} \approx t_{\text{ff,glob}}/8000$. Thus, over the course of Phase II, the inner regions are in a quasi-steady state given the large-scale flow, while the large-scale flow cannot change appreciably, since Phase II lasts for $\lesssim t_{\text{ff,glob}}$.

A feature of our setup is that the CPD has to build up from the infalling gas, since we do not initially put in any structure. The accreting gas naturally accumulates in a CPD whose outer edge is defined by a radial shock and increases over timescales of hundreds of $t_{\text{ff,glob}}$. The thickness (aspect ratio) of the CPD does not vary much while it grows. The formation of the CPD causes a spherical shock to propagate outward through the infalling material, with part of shock at the position of the expanding outer edge of the CPD. Consequently, the flow pattern close to the planet is not quite in steady state. We estimate in Appendix C how much this affects our analysis and find that it should not change our conclusions. Therefore, for simplicity, we will call the Phase II state with a qualitatively constant flow pattern a quasi-steady state and analyze it,

keeping in mind that over much longer timescales, there are likely quantitative fluctuations.

### 3.3. Free Parameters

One can parameterize the degrees of freedom of the problem in different ways. To help bridge simulations and observations, we choose as independent parameters

$$(M_{\star}, a, \Sigma, h, r_{\text{max}}, M_{\text{p}}, R_{\text{p}}), \qquad (7)$$

where, repeating some definitions, $M_{\star}$ is the stellar mass, $a$ is the semimajor axis of the planet, $\Sigma$ is the PPD surface density at $a$ as reduced by the gap opening, $h$ is the PPD aspect ratio at $a$, $M_{\text{p}}$ is the planet mass, and $R_{\text{p}}$ is the physical radius of the planet. Instead of setting $\Sigma$ directly, one could choose a value for the viscosity parameter $\alpha$ (Shakura & Sunyaev 1973) and an unperturbed surface density of the PPD $\Sigma_0$ and, following Kanagawa et al. (2018), let $\Sigma = \Sigma_0/(1+K/25)$, where $K = q^2/(\alpha h^5)$.

To be consistent with the approximation of symmetry around the planet, one should choose $r_{\text{max}}$ to be smaller than the width of the gap, which we do not model. This choice also predicts the reduced surface density to be constant across the gap (Kanagawa et al. 2017), which we assume when setting $\rho(r_{\text{max}})$ (Equation (1)). Also, we introduced the parameter $r_{\text{max}}$ because we do not simulate the whole PPD. However, $r_{\text{max}}$ should be primarily seen not as a numerical parameter but rather as a (simple) way of controlling the incoming angular momentum of the gas. We vary $r_{\text{max}}$ in Section 4.3.

The other characteristic quantities follow from Expression (7): the planet-star mass ratio $q = M_{\text{p}}/M_{\star}$; Hill radius $R_{\text{Hill}} = a(M_{\text{p}}/[3M_{\star}])^{1/3}$; Bondi radius $R_{\text{Bondi}} = aq/h^2$; pressure scale height $H_{\text{P}} = ah$; Keplerian orbital angular frequency of the planet $\Omega_0$ (defined in Equation (3)); and, in particular, $q_{\text{th}}$, the ratio between the planet mass and the disk thermal mass (for short, “thermal mass”; e.g., Korycansky & Papaloizou 1996; Machida & Kokubo 2008; Fung et al. 2019):

$$q_{\text{th}} = q/h^3 = R_{\text{Bondi}}/H_{\text{P}} = 3(R_{\text{Hill}}/H_{\text{P}})^3. \qquad (8)$$

When $R_{\text{p}} \ll \min(R_{\text{Bondi}}, R_{\text{Hill}})$, in the isothermal and inviscid limit, one may expect $q_{\text{th}}$ to be the only parameter controlling the flow in a local region around a planet on a Keplerian orbit (Korycansky & Papaloizou 1996; see, however, Béthune & Rafikov 2019a, 2019b). The radiation transfer introduces a physical scale through the temperature- and density-dependent opacities, but qualitatively, $q_{\text{th}}$ should be key in determining the flow.

One characteristic quantity emerges from our setup: the net mass inflow rate into the Hill sphere $\dot{M}_{\text{Hill, net}}$, that is, what flows in minus what flows out. This in turn is set by more global PPD physics (e.g., Choksi et al. 2023; Nelson et al. 2023). The growth rate of the planet cannot be controlled directly but is at most $\dot{M}_{\text{Hill, net}}$; it is less if some of the large-scale flow instead feeds the CPD. We need to measure $\dot{M}_{\text{Hill, net}}$ from the simulation output because we only set the gradient of the radial velocity at $r_{\text{max}}$ (Section 3.2.2), so that we do not know a priori how much mass will flow in or out as a function of angle. However, the gas at $r_{\text{max}}$ will turn out to be in (inward) freefall at all angles. Then, $\dot{M}_{\text{Hill, net}}$ is maximal and given by $\dot{M}_{\text{Hill, net}} = \int 4\pi r^2 \rho |v_{\text{ff,\infty}}| \sin \theta d\theta$.

<footnote n="6">⁶ By default in Pluto, all cells were considered; we changed this.</footnote>
The flow patterns that we will obtain should not depend sensitively on our choice of $Σ$ and hence $M_{Hill, net}$. This would hold exactly in pure hydrodynamics simulations, but here the optical depth introduces a length scale. However, in practice, this is not an important effect, since the radiative transfer and hence the thermodynamics in the freefall flow do not depend strongly on the density, and even large variations in the Rosseland optical depth do not modify the flow, at least in 1D (Paper I).

### 3.4. Parameter Values Guided by PDS 70 b

We consider parameters that could be appropriate for PDS 70 b (Keppler et al. 2018; Bae et al. 2019; Toci et al. 2020; Wang et al. 2021) without, however, attempting to match the observational properties exactly. We take $M_p = 2 M_J$, guided by the posterior distribution of Wang et al. (2021) and other tentative indications of a low, few-$M_J$ mass (Bae et al. 2019; Stolker et al. 2020; Uyama et al. 2021). A higher value is also possible and is considered in Section 4.4. The surface density is set to $Σ = 0.021 \text{ g cm}^{-2}$ for the gas in the gap of the background PPD at $a = 22 \text{ au}$, coming from $Σ_0 = 2.7 \text{ g cm}^{-2}$ with $α = 10^{-3}$ following<fnref n="7" /> Bae et al. (2019). Bae et al. (2019), or Toci et al. (2020) with their $α ≈ 0.005$, obtain surface densities closer to $Σ ≈ 10^{-3} \text{ g cm}^{-2}$ due to their higher $M_p$. We return to $Σ$ in Section 5.1. Using the expressions of Kanagawa et al. (2017), the bottom of the gap, with a constant surface density, is $Δa_{gap} = 7.1 \text{ au}$ wide, centered on $a = 22 \text{ au}$. Our assumption of a constant $Σ$ over the outer boundary of the simulation domain is thus justified, since $R_{Hill} = 4100 R_J = 2.0 \text{ au}$ is smaller than $Δa_{gap}/2$.

In Table 2, we summarize our choices and the resulting relevant quantities, including the PPD pressure scale height, freefall time, and disk thermal mass of the planet. We are in the high-mass regime with $q_{th} = 7.1 ≫ 1$. This high value of $q_{th}$ is the same as in one of the simulations of Maeda et al. (2022), who, however, used a very different setup (Table 1) and did not study the accretion close to and at the surface of the planet.

For the fiducial case, we let Phase I run for $7.0 × 10^8 \text{ s}$ ($2.0 t_{ff,glob}$) before switching to Phase II. The snapshot used for the analysis was taken at $2.4 × 10^8 \text{ s}$ after the beginning of Phase II, which represents $0.7 t_{ff,glob}$ but thousands of freefall times from $r_{frz} = 10 R_J = 5 R_p$ to $R_p$. As a check, we kept separately running the simulation of Phase I up to $t ∼ 30 t_{ff,glob}$. As expected, the overall flow remained the same, while the CPD grew in size and slightly in thickness. Therefore, the structure of the flow and the CPD in Phase II are representative of a possible steady state. A similar description—several $t_{ff,glob}$ for Phase I and snapshots taken at more than hundreds of freefall times from $r_{frz}$ in Phase II—also applies qualitatively to simulations varying $M_p$ and $r_{max}$ (Sections 4.3 and 4.4).

## 4. Results

Here we present the flow of the gas from the Hill sphere down to the planet and CPD. In Section 4.1, we analyze what fraction of the gas entering the Hill sphere reaches the planet directly and what fraction has sufficient velocity to generate $Hα$. In Section 4.2, we estimate the resulting $Hα$ luminosity and compare it to the assumption that the preshock velocity is $v_{ff}$. In Section 4.3, we assess the effect of the 2.5D

<table><caption>**Table 2**
Chosen and Derived Fiducial Parameters</caption><thead><tr><th>Quantity</th><th>Symbol and Value</th></tr></thead><tbody><tr><td colspan="2"><strong>Chosen Parameters (Expression (7))</strong></td></tr><tr><td>Stellar mass</td><td>$M_* = 0.9 M_☉$</td></tr><tr><td>Semimajor axis</td><td>$a = 22 \text{ au}$</td></tr><tr><td>PPD surface density in gap</td><td>$Σ = 0.021 \text{ g cm}^{-2}$</td></tr><tr><td>PPD aspect ratio at <i>a</i></td><td>$h = 0.067$</td></tr><tr><td>Outer radius of domain</td><td>$r_{max} = 1 \times R_{Hill}$</td></tr><tr><td>Planet mass</td><td>$M_p = 2 M_J$</td></tr><tr><td>Planet radius</td><td>$R_p = 2 R_J$</td></tr><tr><td colspan="2"><strong>Derived Parameters</strong></td></tr><tr><td>Disk thermal mass</td><td>$q_{th} = 3(R_{Hill}/H_p)^3 = 7.1$</td></tr><tr><td>Hill radius</td><td>$R_{Hill} = 4100 R_J = 1.33 H_p$</td></tr><tr><td>Bondi radius</td><td>$R_{Bondi} = 22 \text{ kR}_J = 10.4 \text{ au}$</td></tr><tr><td>Orbital period</td><td>$2\pi/\Omega_0 = 3.4 \times 10^9 \text{ s}$</td></tr><tr><td>Freefall time from $r_{max}$</td><td>$t_{ff,glob} = 3.5 \times 10^8 \text{ s}$</td></tr><tr><td>Mass flux into Hill sphere</td><td>$\dot{M}_{Hill, net} = 6.9 \times 10^{-6} M_J \text{ yr}^{-1}$</td></tr><tr><td>Midplane density in gap</td><td>$\rho_{mid} = 3.8 \times 10^{-16} \text{ g cm}^{-3}$</td></tr><tr><td>Gap 50% full width</td><td>$Δa_{gap} = 7.1 \text{ au}$</td></tr><tr><td>Freefall velocity on planet</td><td>$v_{ff}(R_p) = 59.5 \text{ km s}^{-1}$</td></tr><tr><td colspan="2"><strong>Other Chosen Parameters</strong></td></tr><tr><td>Planet normalized spin</td><td>$f_\omega = \omega/\omega_{crit} = 0.1$</td></tr><tr><td>Maximal dust fraction</td><td>$f_d/g = 0.0001 = 0.01 \times \text{ISM}$</td></tr></tbody></table>

**Note.** To set $Σ$, we chose $Σ_0$ and $α$ from Bae et al. (2019), and the gap width follows from Kanagawa et al. (2017; see text).

approximation by varying the angular momentum of the incoming gas, and in Section 4.4, we consider a higher planet mass.

### 4.1. Gas Flow from the Hill Sphere to the Planet

Figure 1 shows the large- and small-scale structure of the flow and the gas density. At the Hill sphere, the gas is in radial freefall at all latitudes, so that all the gas entering the Hill sphere will accrete onto (that is, become part of) the planet or CPD. For comparison, in 3D, a fraction of the flow would flow back out on perturbed horseshoe orbits (Machida & Kokubo 2008; Lambrechts & Lega 2017; Maeda et al. 2022). Therefore, our inflow rate corresponds to the net inflow in 3D (see also Section 3.3).

Everywhere outside of the planet and CPD (i.e., where the gas is infalling), the radiative flux<fnref n="8" /> is almost purely radial (see Appendix D). Accordingly, the temperature is approximately constant along $θ$ at a given radius. Then, the density stratification leads to a positive pressure gradient: $dP/dθ ∝ dρ/dθ ∝ -dρ/dz > 0$. In turn, the pressure gradient pushes the gas outside of $r ≈ 100 R_J$ poleward. It does so by at most $5°$ compared to a radial trajectory, and inside of $r ≈ 100 R_J$, angular momentum conservation deflects the gas outward. This deviation from a radial trajectory becomes more important closer to the planet. For example, the streamline that started at $θ_{init} ≡ θ(r_{max}) = 20°$ at $R_{Hill}$ joins the CPD at $θ = 73° ≫ θ_{init}$ (Figure 1(c)). The gravitational potential energy of the gas serves to increase all three components of the velocity.

<footnote n="8">⁸ Since we perform radiation hydrodynamical simulations, this automatically includes the interior fluxes from the planet and the CPD, as well as the accretion luminosities (the kinetic energy transformed into radiation at the respective shocks).</footnote>

<footnote n="7">⁷ How this fits with the proposed age for the star of 8–10 Myr (Žerjal et al. 2023) instead of 5 Myr (Müller et al. 2018) should be reassessed.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Density structure (color) and flow pattern (lines) from R<sub>Hill</sub> to R<sub>p</sub> scales. The streamlines start at R<sub>Hill</sub> at θ<sub>init</sub> ≡ θ(r<sub>max</sub>) = 0° 5 and from 5° up to 85° in steps of 5° (panels (a)–(d); thick: multiples of 10°) or also 1° (panels (b) and (d)). Due to angular momentum conservation, the arrow-bearing streamline with θ<sub>init</sub> = 20°, for example, hits the CPD with θ ≫ 20° and not on the planet surface. Only gas within θ<sub>init</sub> ≈ 10° shocks with sufficient velocity to generate Hα (θ<sub>init</sub> = 0° 5, 5°, 10°; black).</figcaption></figure>

The key result seen in Figure 1 is that most streamlines reach the CPD at a large distance (hundreds of Jupiter radii) from the planet. Only a small fraction of the total mass influx from R<sub>Hill</sub> falls in close to the planet. This consequence of angular momentum conservation was seen by Tanigawa et al. (2012) with a very different setup (Table 1). Independent work by Z. Chen & X. Bai (2023, in preparation) also finds this. It has been derived analytically for ballistic (pressure-free) trajectories starting from an outer edge in solid-body rotation in the context of star formation (Ulrich 1976; Mendoza et al. 2009). We have now obtained that most gas reaches the CPD far from the planet when radiation transfer and thermal effects are also included. This conclusion will be seen to hold for other parameter combinations.

Two partial mass influx rates are of interest, in particular for simulations that cannot resolve down to these scales. One is the gas falling directly onto the planetary surface,

$$ \dot{M}_{\text{p, direct}} = 4\pi R_{\text{p}}^2 \int \rho(\theta) v_r(\theta) \sin\theta \, d\theta, \quad (9) $$

integrated at $r = R_p$ from the pole ($\theta = 0$) down to the surface of the CPD, where it connects to the planetary surface (in a so-called “boundary layer”; e.g., Kley 1989; Hertfelder & Kley 2017). This

gives $\dot{M}_{\text{p, direct}} = 3.0 \times 10^{-8} M_J \text{yr}^{-1}$ at $R_p = 2 R_J$, where the CPD is roughly 10° thick. Since $\dot{M}_{\text{Hill, net}} = 6.9 \times 10^{-6} M_J \text{yr}^{-1}$ (Table 2), only $\dot{M}_{\text{p, direct}}/\dot{M}_{\text{Hill, net}} = 0.4\%$ of the gas entering the Hill sphere reaches the planetary surface directly. The expressions of Adams & Batygin (2022) predict a qualitatively similar result, with differences because they assumed a uniform density and solid-body rotation at the Hill sphere and neglected pressure forces (see our Section 4.1), as in Mendoza et al. (2009).

The second partial mass flow rate, $\dot{M}_{\text{Hα}}$, measures the mass influx able to generate Hα, that is, the gas whose preshock velocity $v_0$ exceeds $v_{\text{Hα, crit}}$. To determine $\dot{M}_{\text{Hα}}$, we look at the component of the velocity that is perpendicular to the planetary and CPD surface, shown in Figure 2. For the shock at the surface of the planet, the normal component is the radial velocity $v_r$ because the planet surface is nearly spherical. Even for this small mass of $M_p = 2 M_J$, over the whole free surface of the planet, the gas is fast enough to generate Hα (Equation (2)). Out to at least ≈10 $R_J$ for the CPD in our simulation, the CPD (shock) surface is flat and nearly radial, such that the normal component is nearly equal to the polar velocity $v_θ$ evaluated above the shock surface. For simplicity, we take $v_0 = v_θ$. We identify the largest radius, $r_{30 \text{ km s}^{-1}}$, on the CPD out to which
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Components of the velocity (in kilometers per second) perpendicular to each shock $v_0$: radial $v_r$ around streamlines that will hit the planet surface and polar $v_\theta$ for the ones that will hit the CPD. The Hα emission requires $v_0 > v_{H\alpha, crit} \approx 30$ km s⁻¹ (solid gray streamlines; black dotted otherwise). Streamlines start at $R_{H\alpha}$ in steps of $\Delta\theta = 1^\circ$. Left: fiducial simulation ($M_p = 2 M_J$); right: HigherMass ($M_p = 5 M_J$). The two velocity scales differ, but both have colors only above $v_{H\alpha, crit}$. In both cases, the preshock velocity over the whole free planetary surface is high enough for Hα emission, but only a part of the CPD emits.</figcaption></figure>

$v_0 \ge v_{H\alpha, crit}$ and integrate the mass flux at that radius from the pole down to the surface of the CPD:

$$ \dot{M}_{H\alpha} = 4\pi r_{30 \text{ km s}^{-1}}^2 \int \rho(\theta) v_r(\theta) \sin\theta \, d\theta. \quad (10) $$

Due to time independence, this is equivalent to an integral over the shock surfaces. We find a maximum radius $r_{30 \text{ km s}^{-1}} = 3.3 R_J$ and a CPD height of $14^\circ$ there. This yields $\dot{M}_{H\alpha} = 4.7 \times 10^{-8} M_J \text{ yr}^{-1}$, or $\dot{M}_{H\alpha}/\dot{M}_{H\alpha, net} = 0.7\%$.

Both $\dot{M}_{p, direct}$ and $\dot{M}_{H\alpha}$ are small fractions of $\dot{M}_{H\alpha, net}$. Correspondingly, they originate from a narrow polar region in which the specific angular momentum of the gas $j_z = Rv_\phi \propto r^2 \sin^2\theta$ (Equation (3)) is low. Indeed, tracing the streamlines that define $\dot{M}_{p, direct}$ and $\dot{M}_{H\alpha}$ back to $r_{max}$, we find starting angles of $\theta_{init} \approx 7^\circ$ and $9^\circ$, as seen in Figure 2.

Comparing $\dot{M}_{H\alpha}$ and $\dot{M}_{p, direct}$, we see that $(1-\dot{M}_{p, direct}/\dot{M}_{H\alpha}) \approx 60\%$ of the total Hα-generating gas falls on the CPD surface and not on the planet. This assumes that all of $\dot{M}_{p, direct}$ produces Hα, which holds. However, the local Hα flux $F_{H\alpha}$ depends strongly on $v_0$ (crudely, $F_{H\alpha} \sim v_0^3$; Aoyama et al. 2018), so that it is not clear a priori whether the CPD or the planetary surface dominates the total emission. We look at this in more detail in the next section.

## 4.2. Approximate Hα Emission

We estimate the observable Hα luminosity from the planet and CPD surface shocks. Detailed 2D radiation transport of the generated Hα is beyond the scope of this paper. However, based on Marleau et al. (2022) and for $f_{d/g} = 10^{-4}$, we expect the incoming gas and dust to be very optically thin to Hα photons<fnref n="9" /> for a mass influx rate of $\dot{M}_{H\alpha, net} \sim 10^{-5} M_J \text{ yr}^{-1}$ or even much higher. Also, given that $q_{th} \gtrsim 1$ planets open gaps, extinction by the PPD is possibly negligible. Therefore, summing the local Hα production along the planetary and CPD surfaces (the radiative source terms) gives a reasonable estimate of the luminosity leaving the system.

We display in Figure 3 the Hα flux per emitting area $F_{H\alpha}$, which depends only on the local preshock density $\rho_0$ and preshock velocity $v_0$ above either shock. We use the data of Aoyama et al. (2018) for $F_{H\alpha}(\rho_0, v_0)$. For the surface shock, the flux is plotted against the distance from the pole, and for the CPD shock, it is plotted outward from the surface of the planet, where the CPD begins. Instead of the linear distance, we use

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Flux from the accretion shocks on the planet (solid red curve and shaded area) and CPD (blue) surface as a function of the cumulative area along the planetary surface starting at the pole or the CPD surface starting at the planet at $R_p = 2 R_J$ (Equation 11(a)). With this x-axis, the area under each curve is proportional to its contribution to the luminosity. The spike comes from the high-velocity surface flow but barely contributes to the total flux, contrary to what Takasao et al. (2021) found (see Section 5.2). Dashed curves show the cumulative integral of each contribution (right axis).</figcaption></figure>

the respective cumulative areas (including both hemispheres),

$$ A_{\text{plnt surf}}(\theta) = \int_{0}^{\theta} 4\pi R_{p}^{2} \sin\theta' d\theta', \quad (11a) $$

$$ A_{\text{CPD surf}}(r) = \int_{R_p}^{r} 4\pi r' \sin\theta \, dr', \quad (11b) $$

where we have assumed a $\theta = \text{constant}$ CPD surface in Equation 11(b), which holds approximately for the region whose Hα emission dominates (see Figure 2). For the actual analysis, we look for the temperature peak (the Zel'dovich spike) for each $r = \text{constant}$ ring in the $r-\theta$ plane and use the cell above it to define the surface. With Equation 11(a) to measure distance, the ratio of the areas under the curves gives the relative contribution of each shock. Emission comes from the exposed planetary surface from the pole down to the surface of the CPD at $\theta_{max}$, which corresponds to a filling factor $f_{fill} = 1 - \cos\theta_{max} = 0.83$, and from the CPD from the planetary radius out to $r \approx r_{30 \text{ km s}^{-1}} \approx 3.3 R_J$ (as found in Section 4.1). The region without emission is labeled in Figure 2.

The integrated luminosities (dashed lines in Figure 3) are $L_{H\alpha} = 3.2 \times 10^{-8} L_\odot$ for the planetary surface shock and

<footnote n="9">⁹ However, the Planck mean opacity is high enough for the gas and radiation temperatures to be equal, except in the Zel'dovich spikes (Appendix D).</footnote>
$L_{H\alpha} = 0.45 \times 10^{-8} L_{\odot}$ for the contribution by the CPD. Approximately, ignoring the angular dependence of the radiation by summing the two terms, an observer looking at the system would see an $H\alpha$ luminosity of $L_{H\alpha} \approx 3.7 \times 10^{-8} L_{\odot}$. We compare this to the observations of PDS 70 b in Section 5.1. The supersonic surface flow on the CPD, discussed in Appendix B, causes a local spike in the emission (shaded region in Figure 3). However, this thin layer contributes negligibly to the integrated emission. This contrasts strongly with the results of Takasao et al. (2021), to which we return in Section 5.2.

Figure 3 shows that the planetary surface shock largely dominates the $H\alpha$ emission. The preshock densities, $\rho_0 \sim 10^{-13} \text{ g cm}^{-3}$, which depend only weakly on position, are similar to within 0.1 dex between both shocks, and the emitting areas are similar (near $45 R_p^2$). However, the difference in the preshock velocities is much more consequential. The velocities are different in part because the emitting region of the CPD is at a slightly greater $r$ (Figure 2; $r \approx 3$ and (3–10) $R_J$) for the two simulations) than the planet surface ($R_p = 2 R_J$). The other, and more important, factor is that for the CPD, the shock velocity (i.e., the component orthogonal to the shock surface) is the polar velocity $v_\theta$, and this is even smaller than the local radial velocity $v_r$.

We can compare the $H\alpha$ luminosity from the free planetary surface to the luminosity expected from purely radial accretion at $v_{ff}$ for the same $\dot{M}_{\text{p, direct}}$ and $f_{\text{fill}}$. Explicitly, the combination ($\dot{M}_{\text{p, direct}}$, $M_{\text{p}}$, $R_{\text{p}}$, $f_{\text{fill}}$) that we have here implies an average preshock number density for the gas right above the shock $n_0 = X \dot{M}_{\text{p, direct}} / (4\pi R_p^2 f_{\text{fill}} v_{ff} m_{\text{H}}) = 6.3 \times 10^{10} \text{ cm}^{-3}$ (Equation (A9) of Aoyama et al. 2020) and thus $L_{H\alpha}' = 4\pi R_p^2 f_{\text{fill}} \times F_{H\alpha}(n_0, v_{ff}) = 4.7 \times 10^{-8} L_{\odot}$, which is higher by 50% than what we found. The difference is due mainly to the strong dependence of $F_{H\alpha}$ on $v_0$ with roughly $F_{H\alpha} \propto v_0^3$ (Aoyama et al. 2018). Indeed, in our simulation, the radial velocity at the pole is equal to the freefall value $v_{ff} = 59.5 \text{ km s}^{-1}$, but at the equator, it is lower by about 30%. This is because the gas gains more velocity in $\theta$ and $\phi$ further away from the pole, with the three components summing up to $v_r^2 + v_\theta^2 + v_\phi^2 = v_{ff}^2$ everywhere in the freefalling region by conservation of energy. In other words, centrifugal forces due to angular momentum conservation are slowing down the infalling gas at low latitudes. This reduction of the preshock velocity leads to less emission ($\sim \rho_0 v_0^3$) for the same total mass flow rate ($\sim \rho_0 v_0$) compared to the simple assumption of radial freefall.

### 4.3. Varying the Incoming Angular Momentum

The specific angular momentum of the gas entering the domain depends on the choice of the outer radius $r_{\text{max}}$ (Equation (3)). As argued in Section 3.1, setting $r_{\text{max}} = R_{\text{Hill}}$ is a natural choice for simulations assuming axisymmetry around the planet, but it remains approximate. The incoming angular momentum sets what fraction of the gas can reach the planet, and in general, it might determine whether an outflow near the midplane occurs or not. Here we have infall at all angles, but an outflow does occur in the azimuthal average of the 3D simulations (e.g., Schulik et al. 2020).

Therefore, we performed two additional simulations with the same parameters except for $r_{\text{max}} = 0.7$ and 1.3 $R_{\text{Hill}}$, named LowAngMom and HighAngMom, to give the accreting gas less or more angular momentum, respectively. We chose planet

radii near $R_p = 2 R_J$ by keeping $r_{\text{min}} = 1.9 R_J$, with the exact values set by how much mass is accreted and how it cools until a quasi-steady state is established as described in Section 3.2.4. The radii turn out to be, respectively, $R_p = 2.15$ and $1.95 R_J$, which is similar enough for our purposes, especially since the radius does not directly affect the infall of matter (see also Section 5.3). Since $R_{\text{Hill}} = 1.33 H_p$ (Table 2), LowAngMom has $r_{\text{max}} \approx H_p$. We took $N_\theta = 51$ instead of $N_\theta = 181$, since it does not influence the accretion flow. We also find that for these simulations, the radial velocity of the gas at $r_{\text{max}}$ quickly reaches and remains at the freefall velocity for all angles. Because we keep the density at $r_{\text{max}}$ the same, the mass influxes at the Hill sphere are somewhat lower and higher, respectively, with $\dot{M}_{\text{Hill, net}} = 5 \times 10^{-6} M_J \text{ yr}^{-1}$ for LowAngMom and $\dot{M}_{\text{Hill, net}} = 9 \times 10^{-6} M_J \text{ yr}^{-1}$ for HighAngMom, instead of $\dot{M}_{\text{Hill, net}} = 7 \times 10^{-6} M_J \text{ yr}^{-1}$ for the fiducial run. As summarized in Section 3.3, judging from Paper I, this difference in density will not influence the flow pattern.

Figure 4 compares the flow of the gas on scales of $r \sim 50–150 R_J \sim 0.01–0.03 R_{\text{Hill}}$ for $r_{\text{max}}/R_{\text{Hill}} = 0.7, 1.0$, and 1.3. As expected, with larger $r_{\text{max}}$, only gas from a smaller cone around the pole reaches the planetary surface directly (highlighted in pink). The starting angle of the last streamline hitting the planet surface is $\theta_{\text{init}} = 15^\circ$8 for LowAngMom and $\theta_{\text{init}} = 3^\circ$9 for HighAngMom, which bracket the corresponding $\theta_{\text{init}} = 7^\circ$1 for the fiducial case. The same applies to the mass fluxes relative to the respective $\dot{M}_{\text{Hill, net}}$, which are $\dot{M}_{\text{p, direct}}/\dot{M}_{\text{Hill, net}} = 2.4\%$ for the simulation LowAngMom and $\dot{M}_{\text{p, direct}}/\dot{M}_{\text{Hill, net}} = 0.052\%$ for HighAngMom; the fiducial case had $\dot{M}_{\text{p, direct}}/\dot{M}_{\text{Hill, net}} = 0.7\%$.

Thus, at fixed $\dot{M}_{\text{Hill, net}}$, $\dot{M}_{\text{p, direct}}$ depends sensitively on $r_{\text{max}}$, but, overall, the fraction of $\dot{M}_{\text{Hill, net}}$ that shocks on the planetary surface is at most of order of 1%. The dependence of $\dot{M}_{\text{Hill, net}}$ on $r_{\text{max}}/R_{\text{Hill}}$ will be similar. We discuss this further in Section 5.3. For reference, we obtain $L_{H\alpha} = 1.0 \times 10^{-7}$ and $7.3 \times 10^{-9} L_{\odot}$ for LowAngMom and HighAngMom, respectively, again bracketing the fiducial case with its $L_{H\alpha} = 3.7 \times 10^{-8} L_{\odot}$.

### 4.4. Varying the Planet Mass

Short of doing a full exploration of the whole parameter space, we consider a higher $M_{\text{p}}$. The planetary mass is an important parameter that controls the dynamics of the gas. We wish to see whether here too only a fraction of the gas entering the Hill sphere falls directly onto the planet and whether the planetary surface still dominates the $H\alpha$ emission relative to the CPD, as it does for $M_{\text{p}} = 2 M_J$ (Figure 3). A priori, especially the latter could change for higher masses because $v_0$ scales to zeroth order with $v_{ff} \propto r^{-1/2}$, which will lead to a larger $H\alpha$-emitting area on the CPD.

We therefore simulate an accreting planet as in the fiducial run but with $M_{\text{p}} = 5 M_J$ and name this HigherMass. Bae et al. (2019) used this value for PDS 70 b. We again set $r_{\text{max}}$ equal to $R_{\text{Hill}}(M_{\text{p}})$ and use $N_\theta = 51$. The other parameters are the same, leading to a large $q_{\text{th}} = 18$. Figure 5 shows that the flow pattern for the $5 M_J$ simulation is similar to the one of the $2 M_J$ simulation. Even at this lower resolution in $\theta$, a thin ($\sim 1^\circ$) fast inward surface flow is seen again. Repeating the analysis above, we find that for $5 M_J$, only $\dot{M}_{\text{p, direct}}/\dot{M}_{\text{Hill, net}} = 1.2\%$ of $\dot{M}_{\text{Hill, net}} = 1.4 \times 10^{-5} M_J \text{ yr}^{-1}$ reaches the planet directly, and that $\dot{M}_{H\alpha}/\dot{M}_{\text{Hill, net}} = 7\%$ shocks with $v \gtrsim v_{H\alpha}$, crit. For
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Fraction of the total inflow rate reaching the planet surface for different amounts of angular momentum at the Hill sphere, probed by varying the position of the outer radius of the simulation domain: $r_{max}/R_{Hill} = 0.7, 1.0$ (fiducial run), and 1.3 (left, middle, and right panel, respectively). In all cases, the gas flows inward for all angles from $r_{max}$ down to the planet or CPD. The gray scale shows the density (logarithmic). Each simulation extends to $r_{max} = 4100 R_J$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of the flow in the $M_p = 5 M_J$ simulation (black) to the fiducial 2 $M_J$ run (red). Streamlines are for $\theta_{\text{init}} = 0^\circ, 5^\circ, \dots, 90^\circ$ (left) or $\theta_{\text{init}} = 0^\circ, 1^\circ, 2^\circ, \dots$ (right). The axes, different for each panel, are scaled by the respective Hill radius of each simulation. Dark blue dots show the cell centers in the 5 $M_J$ run.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Same as Figure 3 but for 5 $M_J$. The total is $L_{H\alpha} = 1.4 \times 10^{-7} L_{\odot}$. The planetary surface is twice as bright as the CPD.</figcaption></figure>

$M_p = 2 M_J$, we recall that we had smaller fractions of 0.4% and 0.7%, respectively.

In Figure 6, we show the $H\alpha$ emission from the planetary surface and the CPD, as in Figure 3. Again, the inward flow below the CPD surface does not generate an overall important $H\alpha$ flux; there is a local spike, but its relative contribution is negligible. The planetary surface generates $L_{H\alpha} = 1.2 \times 10^{-7} L_{\odot}$, and the CPD surface generates $L_{H\alpha} = 0.54 \times 10^{-7} L_{\odot}$. Summing the two terms again in lieu of detailed radiation transport yields a total of $L_{H\alpha} \approx 1.7 \times 10^{-7} L_{\odot}$. Thus, the CPD surface emits about 30% of the total flux, up from 15% in the 2 $M_J$ case. This relative increase is because $v_0 > v_{H\alpha, \text{crit}}$ out to a cylindrical radius of $R \approx 9 R_J$ instead of $R \approx 3 R_J$ in the fiducial case (see the colored regions in both panels of Figure 2).

## 5. Discussion

Our main results are that (a) only a small fraction of the net mass flux into the Hill sphere falls directly onto the planet, (b) only a slightly larger fraction produces any $H\alpha$ (Figure 1 and 2), and (c) the emitted $H\alpha$ comes from both the planetary and CPD surfaces and not from the fast flow beneath the CPD surface (Figure 2). Our simulations were conducted in 2.5D, but Tanigawa et al. (2012) obtained qualitatively the same Hill sphere flow structure in their larger-scale isothermal 3D simulations. This lends support to our 2.5D approach and suggests that these findings are robust. The advantage of 2.5D is that it makes it computationally much more accessible not to smooth the gravitational potential while including radiation transport. This allowed us to simulate down to subplanet scales, crucial for calculating $L_{H\alpha}$, since the highest-velocity regions strongly dominate the emission.

In Section 5.1, we look at PDS 70 b. In Section 5.2, we compare our results with other predictions of $H\alpha$ emission from planets accreting other than by magnetospheric accretion (for the latter, see the discussion in Aoyama et al. 2021). Finally, in Section 5.3, we comment on a few aspects of our models.

### 5.1. Comparison with PDS 70 b

As a check, we compare the $L_{H\alpha}$ estimated from the simulations with the observational data for PDS 70 b. This planet had motivated our parameter choices (Table 2).

Assuming that the $H\alpha$ photons are leaving the system isotropically, its measured luminosity is $L_{H\alpha} = 7 \times 10^{-7} L_{\odot}$ (Zhou et al. 2021; Sanghi et al. 2022). At such luminosities, absorption within the system is likely unimportant for a very wide range of dust opacities (Marleau et al. 2022), and absorption by the PPD is more likely to be low given that the planet is found in a gap. Therefore, a direct comparison is meaningful.

For the different runs, we obtained $L_{H\alpha} \sim 10^{-8}–10^{-7} L_{\odot}$, which is 10–100 times smaller than the observationally derived value. This is in fact satisfactory given that we took nominal model parameters (Table 2) from the literature without efforts to match the $L_{H\alpha}$. Reducing the incoming angular momentum or radius of the planet or using $M_p = 12 M_J$ as in Dong et al. (2021) would make it easy to raise our $L_{H\alpha}$ closer to the measured value. Beyond this, magnetospheric accretion columns, if present, could also be contributing to the flux.

### 5.2. Comparison to Other Predictions of $H\alpha$

Only a few studies so far predict the $H\alpha$ emission of forming planets. Thanathibodee et al. (2019) applied magnetospheric accretion radiation transfer models for stars to the planetary regime, but Szulágyi & Ercolano (2020) were the first to present hydrogen-line luminosities based on 3D radiation hydrodynamics simulations. However, their smoothing of the gravitational potential (Table 1) makes their results challenging to interpret, as Aoyama et al. (2020, 2021) discussed. Therefore, we restrict our comparison here to the work of Takasao et al. (2021), who also set $\epsilon_{grav} = 0$.

We follow a similar approach to Takasao et al. (2021) to calculate the $H\alpha$ emission from the simulation data by integrating the local $H\alpha$ flux predicted by Aoyama et al. (2018) as a function of the preshock velocity and density (Section 4.2). However, Takasao et al. (2021) found that nearly 95% of the $H\alpha$ emitted from the planetary surface comes from an approximately 10°–15° thick surface layer of the CPD where the flow hits the planetary surface (see their schematic Figure 13). This differs significantly from our results, in which the sub-CPD surface flow is very thin (~1°) and contributes negligibly to the total luminosity (Figures 3 and 6).

To understand the difference, we compare in Figure 7 the radial velocity and density at $r = R_p$ in the work of Takasao et al. (2021; top panel) and our work (middle panel). Three similar zones are present in both works: the free surface of the planet from the pole down to (geometrically, or “up to” in angle) $\theta = 67^\circ.5$ in their case or $\theta = 80^\circ$ in ours; the fast (sub-) CPD surface flow, which is quite thick in their case and thin in ours (colored regions in Figure 7); and the CPD connecting to the planet surface at $\theta \gtrsim 83^\circ$ for both. In both simulations, the radial velocity above the planet is equal to $v_{ff,\infty}$ at the pole (not shown) and decreases with $\theta$, as discussed in Section 4.2. Also, in both cases, the density increases quickly in the CPD zone, as expected for (approximately) isothermal structures (see Figure 10) in hydrostatic equilibrium.

However, there is a large difference in the layer located between the freefalling gas and the CPD. Takasao et al. (2021) obtained a zone roughly 10°–15° thick in which the radial velocity decreases only slowly by half, which they called a “postshock, converged accretion flow.” However, we find a zone that is only a few degrees thin and in which the radial velocity decreases quickly as a function of distance below the surface. In our simulations, there is also a visible convergence
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Structures in the $\gamma = 1.01$ hydrodynamics simulation of Takasao et al. (2021) compared to ours (with FLD). Top and middle panels: radial velocity immediately above the planet surface, normalized to the respective $v_{\text{If}}$ (Equation (2)) and density (normalized to the maximum, i.e., midplane density) in their (top panel; blue curves) and our (middle panel; pink curves) simulations. Bottom panel: local $H\alpha$ flux times $\sin\theta$, normalized to the respective maxima. This way, the relative areas under the curves are proportional to their contribution to the total flux.</figcaption></figure>

of the postshock accretion flow, seen as the approximately constant-$\theta$ segments of the sub-CPD surface streamlines in Figure 5. However, this converged flow has a smaller $|v_r|$, is thin, and involves relatively little mass. Since the flow is the same while more than doubling the mass ($M_p = 2$ versus $5 M_J$), it most likely does not matter that Takasao et al. (2021) set an even higher mass, $M_p = 12 M_J$.

This qualitative difference in the flows at the CPD surface must come instead from the different thermodynamics that Takasao et al. (2021) assumed, namely, no radiative transfer but an adiabatic EOS with $\gamma$ close to but above unity. Even though $\gamma$ is close to the isothermal value of unity, their adiabatic EOS leads to a thick and hot postshock region surrounding the CPD, as their Figure 6 shows. Consequently, the radial velocity remains high after the shock in $\theta$, so that the gas hits the planet surface quickly; the gas is subsonic, but the Mach number is large, so that the absolute velocity is high. On the contrary, with radiative transfer, the gas cools quickly and the density increases much more across the shock. By mass conservation, the postshock radial velocity is correspondingly smaller, which significantly decreases the amount of emission.

Figure 7(c) shows the local $H\alpha$ line emission at $r = R_p = 2 R_J$ as a function of angle from the pole. We plot $F_{H\alpha} \times \sin(\theta)$, as in Equation (13) of Takasao et al. (2021), and normalize the curves from Takasao et al. (2021) and our simulation independently to their respective maximum. This way, the areas under the curves are proportional to the contribution of each region to the total flux from a model. This shows very clearly that in the simulation of Takasao et al. (2021), only the sub-CPD surface flow, where it hits the planets, generates appreciable amounts of $H\alpha$, whereas in our

case, that zone is negligible for the integral (seen also in Figures 3 and 6).

In our finding that only a small fraction of the large-scale flow falls directly onto the planet, however, we agree qualitatively with Takasao et al. (2021). Excluding the radial flow below the CPD surface, the accretion rate directly onto the planetary surface in their case<fnref n="10" /> is $\dot{M}_{\text{p, direct}} \approx 3 \times 10^{-10} M_J \text{ yr}^{-1}$, which is 0.2% of their net mass influx rate $\dot{M}_{\text{Hill, net}} = 2 \times 10^{-7} M_J \text{ yr}^{-1}$. This is smaller than but similar to our fractions for 2 and 5 $M_J$ (0.4% and 1.2%), with the difference likely coming from their choice of a purely vertical mass flow at their $r_{\text{max}} = 100 R_p$.

We conclude that including radiation transfer in hydrodynamical simulations is important for accurate predictions of $H\alpha$ emission because of their sensitivity to the velocity structure. The flow in the large-scale, supersonic region can likely be well captured by isothermal simulations, but the postshock behavior of the gas depends on the thermodynamics. Smoothing-free 3D simulations in the high-mass (high-$q_{\text{th}}$), low-$R_p/R_{\text{Hill}}$ regime would be a worthwhile complement to the existing work (Table 1).

### 5.3. Further Aspects within and beyond Our Model

We comment on a few aspects within or beyond our model.

**Regions traced by the $H\alpha$.** The results of Section 4.4 suggest that as the planet mass increases, the contribution of the CPD to the total $H\alpha$ becomes increasingly important. However, the two terms remain of the same order of magnitude, and if the CPD contribution is thicker than in Figure 5, the reduced shock velocity would lead to a smaller contribution. Modeling of the line shapes should therefore take both components into account.

**Varying the planetary radius.** We can do this approximately without additional simulations by measuring the $\dot{M}$ according to Equations (9) or (10) at a different $r = R_p'$ and similarly beginning at $R_p'$ the outward integration of the $H\alpha$ emission along the CPD surface. To first order, the choice of $R_p$ will not affect the supersonic flow. In a similar approach, Takasao et al. (2021) set an open boundary at their $r_{\text{min}}$ and used the density and velocity there to calculate the $H\alpha$ emission that would come from a shock at that position. Doing this, we find roughly $L_{H\alpha} \propto 1/R_p$, which can be used to approximately scale the results of one simulation to other $R_p$ values.

**Choice of $r_{\text{max}}$.** We have varied $r_{\text{max}}$ by 30%, with the case $r_{\text{max}} = 1.3 R_{\text{Hill}}$ corresponding to $r_{\text{max}} \approx H_p$. It would be surprising if 3D simulations corresponded effectively to a much larger $r_{\text{max}}$, but the effective $r_{\text{max}}$ could conceivably be smaller. Then, a larger fraction of $\dot{M}_{\text{Hill, net}}$ would reach the planet directly and emit $H\alpha$.

**Models of 1D planet structure.** Global formation models use $\dot{M}_{\text{Hill, net}}$ to set the ram pressure at the surface of the planet when calculating its radius and luminosity (Mordasini et al. 2012b). However, this is an incorrect assumption, since only the much smaller rate $\dot{M}_{\text{p, direct}} \ll \dot{M}_{\text{Hill, net}}$ will set the pressure on the surface of the planet, which could affect its postformation luminosity (e.g., Mordasini 2013; Berardo et al. 2017). In our simulations so far, the (radial) ram

<footnote n="10">¹⁰ We use $\dot{M}_{\text{p, direct}} \approx (1 - \cos\theta_{\text{max}})4\pi R_p^2 \rho_0 v_0$ and read off the values from their Figure 8. Similarly, $\dot{M}_{\text{Hill, net}}$ is close to their maximal $\dot{M}_{\text{Hill, net}}$, which is set by the boundary conditions at their $r_{\text{max}}$ (see their Figure 2(b)).</footnote>
pressure $P_{\text{ram}}(\theta) = \rho(\theta)v_r(\theta)^2$ turns out to be almost constant with polar angle (not shown), so that the reduced ram pressure could be easily included in 1D planet models. However, an appropriate treatment of the boundary layer with its $f_\omega$-dependent transfer of mass, angular momentum, and energy would be needed (e.g., Dong et al. 2021).

*Other hydrogen lines.* Other hydrogen lines such as H$\beta$, Pa$\beta$, or Br$\gamma$ have been observed at a few planetary-mass objects such as Delorme 1(AB)b (Eriksson et al. 2020; Betti et al. 2022a, 2022b). These lines also require a similar minimum shock velocity $v_{\text{crit}} \approx v_{\text{H}\alpha, \text{crit}} = 30 \text{ km s}^{-1}$ to be emitted, since their excitation energies are similar (Aoyama et al. 2018). Therefore, our analysis could have applied to the other lines as well.

*Magnetospheric accretion.* If it proceeds as for young stars (e.g., Romanova et al. 2002; Hartmann et al. 2016), magnetospheric accretion is an interesting mechanism that could let the gas slide ballistically along the magnetic field lines connecting the inner edge of the CPD and the planet surface (Lovelace et al. 2011). This would lead to a shock at the planet surface at almost freefall velocity and thus to H$\alpha$ emission. This would be in addition to what the direct infall $\dot{M}_{\text{p, direct}}$ generates, contrasting with the stellar case, in which $\dot{M}_{\text{p, direct}}$ is essentially zero. In fact, magnetospheric accretion would let almost the same amount of H$\alpha$ be generated as in the 1D spherically symmetric classical picture (e.g., Bodenheimer et al. 2000), since most accreting gas would ultimately reach the planet at (nearly) freefall velocity.<fnref n="11" />

Whether magnetospheric accretion from the CPD onto the planet actually takes place or not is not yet clear. It requires a few conditions to be met: (i) the CPD must be an accretion and not a decretion disk, (ii) the magnetic field of a young planet needs to be able to disrupt the CPD, and (iii) the gas must be sufficiently ionized to couple to the magnetic field (e.g., Keith & Wardle 2014; Hasegawa et al. 2021). In the picture painted by Batygin (2018), in which gas falls toward the pole and a decretion disk, CPD disruption would not be needed, and the apex of the magnetic field lines would increase the effective area of the planetary surface intercepting the flow. So far, interesting but only tentative scaling arguments support the main hypothesis of a sufficiently strong magnetic field (Christensen et al. 2009; Katarzyński et al. 2016), requiring further studies for a robust assessment. Further motivation might come from tentative observational evidence for magnetospheric accretion in the somewhat older, essentially isolated object Delorme 1(AB)b (Ringqvist et al. 2023).

## 6. Summary and Conclusions

We studied the gas flow from the Hill radius down to the surface of a forming super-Jupiter planet able to generate hydrogen lines such as H$\alpha$. We performed axisymmetric, 2D radiation hydrodynamical simulations in a vertical frame centered on the planet and following it on its orbit around the star. These simulations connect to global disk simulations through the net mass inflow into the domain $\dot{M}_{\text{Hill, net}}$ and the angular momentum of the gas, both of which are input parameters here. We argued that the flow structure should depend only a little on $\dot{M}_{\text{Hill, net}}$. Therefore, this should also

apply to the partial accretion rates or the relative contributions to line emission by the planetary and CPD surface.

Two important features compared to previous work are that we (i) included radiation transfer, with tabulated dust and gas opacities, to model correctly the thermal effects that can influence the flow, especially below the CPD surface shock, and (ii) did not smooth the gravitational potential and used a high spatial resolution close to the planetary surface (~0.01 $R_J$). Whereas previous work with a nonzero smoothing length (e.g., Tanigawa et al. 2012) was concerned with the accretion of mass and angular momentum onto the CPD, we focus on the planet surface and the innermost regions of the CPD close to it.

We confirmed that most of the mass flux flowing toward the CPD and the planet forms an accretion shock on the surface of the CPD (Figure 2). Only a very small fraction, of the order of 1%, reaches the planet surface directly, and the fraction shocking at sufficiently high velocity ($v_0 > v_{\text{H}\alpha, \text{crit}} = 30 \text{ km s}^{-1}$) to generate hydrogen lines such as H$\alpha$ is similarly small (Figure 4). We found that these results are robust to variations in the planetary mass and the angular momentum of the incoming gas. The large-scale flow pattern agrees qualitatively with 3D isothermal simulations with a smoothed gravitational potential (Tanigawa et al. 2012; Fung et al. 2019), lending support to our approach.

For all simulations, we estimated the H$\alpha$ emission through the nonequilibrium shock models of Aoyama et al. (2018). Our inclusion of radiative transfer keeps the fast flow beneath the CPD shock surface thin, so that only the free surfaces of the planet and CPD appreciably emit shock tracers. This contrasts with the results of hydrodynamics-only simulations (Section 5.2), showing the importance of including radiation transfer while not smoothing the gravitational potential.

In summary, we have studied one aspect of determining how many planets can be detected at accretion tracers such as H$\alpha$: what parts of the flow can generate accretion-line emission. However, we have not addressed the relation between this H$\alpha$-generating mass flux and the growth rate of the planet. This is a different question, beyond the scope of our work, and involves studying the timescale for mass transport in the CPD. In one limit, only what falls directly onto the planet will let it grow at a given time, but in the other, the CPD would be able to appreciably feed the planet (Adams & Batygin 2022). Dedicated simulations are required.

It is, moreover, a separate issue as to whether the statistics of known accreting planets match expectations given our current understanding of planet formation and the empirical demographics of directly imaged planets (e.g., Nielsen et al. 2019; Vigan et al. 2021). Both the migration and formation timescales influence this, as well as the non-Gaussianity in the residuals in high-contrast images (Marois et al. 2008; see applefy by Bonse et al. 2023). A careful statistical treatment would be welcome (R. Dong et al. 2023, in preparation), as would more detections—for which there is hope, thanks to instrumental progress such as VIS-X (Haffert et al. 2021), KPIC (Delorme et al. 2021), and RISTRETTO (Chazelas et al. 2020), to name a few.

## Acknowledgments

We dedicate this paper to the memory of Willy Kley (Universität Tübingen, †2021 December 21), with whom several of us had the privilege of collaborating and who

<footnote n="11">¹¹ Half of the potential energy is dissipated in the CPD if it extends to the surface of the planet (e.g., Pringle 1981; Hartmann et al. 1997).</footnote>
provided useful and interesting advice at an early stage of this
project. His kindness and expertise are deeply missed. We are
indebted to Yuhiko Aoyama for his model data and for many
answers to many questions. We warmly thank Shinsuke
Takasao, Kazuhiro Kanagawa, and Matthew Bate for helpful
discussions. G.-D.M. is particularly grateful to Thomas
Henning for his vision, support, and encouragement to begin
investigating the accretion shock in planet formation. We thank
the referee for a report that helped us significantly clarify the
structure and content of this paper. This research was supported
by the Munich Institute for Astro-, Particle and BioPhysics
(MIAPbP), which is funded by the Deutsche Forschungsge-
meinschaft (DFG, German Research Foundation) under
Germany's Excellence Strategy (EXC-2094—390783311).
G.-D.M. and R.K. acknowledge the support of the DFG
priority program SPP 1992 "Exploring the Diversity of
Extrasolar Planets" (MA 9185/1, KU 2849/7, and KU 2849/
10). G.-D.M. and C.M. also acknowledge the support from
the Swiss National Science Foundation under grant
200021_204847 "PlanetsInTime." R.K. acknowledges finan-
cial support via the Heisenberg Research Grant funded by the
DFG under grant KU 2849/9. W.B. acknowledges funding by
the DFG under grant KL 650/31-1. Parts of this work have
been carried out within the framework of the NCCR PlanetS
supported by the Swiss National Science Foundation. This
research has made use of NASA's Astrophysics Data System
Bibliographic Services. We gratefully acknowledge the use of
Ankit Rohatgi's WebPlotDigitizer. All figures were
produced using gnuplot with the terminal epslatex and
the font package fouriernc.

**Appendix A**
**Radial Gridding**

The radial gridding is made of three parts and shown in
Figure 8. The inner section at $r \in [r_{\min}, r_{\min} + L_u]$ has 32
uniformly spaced cells $\Delta r_u = 0.001 R_J$ long (hence
$L_u = 0.032 R_J$); the outer section at $r \in [r_{\min} + L_u +
L_s, r_{\max}]$, with $L_s = 0.5 R_J$, is logarithmically stretched with
76 cells per decade in radius; and the transition section at
$r \in [r_{\min} + L_u, r_{\min} + L_u + L_s]$ has geometrically stretched
cells (Mignone et al. 2007) chosen to have a smooth increase
in cell size between $\Delta r_u$ and the first cell size in the logarithmic
part. We take 30 cells for the middle section. This gives 307
zones in total. We have tested that the results do not change
appreciably when using a higher resolution for the different
parts of the grid. As in our 1D simulations (Papers I and II), a
lower resolution in the inner uniform part would lead to
artificially high luminosities in the settling zone below the
shock, with a rapid increase in time. However, for some
simulations, we were able to increase the cell size in the inner
part to $\Delta r_u = 0.002 R_J$ (adjusting the stretched transition region

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Cell size of the radial grid, from $r_{\min} = 1.9 R_J$ to $r_{\max} = R_{\text{Hill}} = 4100 R_J$ in the fiducial case, with its uniform (blue), geometrically stretched (peach), and logarithmic (red) segments. (a) Global view. Only the central position of every fifth cell is shown (vertical lines). (b) Zoom-in near $r_{\min}$. Every cell center is shown.</figcaption></figure>

to have a smooth change in cell size) and still obtain a correct-
looking solution.

**Appendix B**
**Dependence of the CPD Surface Flow on the Resolution in**
**the Polar Direction**

Beneath the shock on the CPD surface, in which the polar
component of the velocity goes from super- to subsonic, there
is a thin layer of ≈1° in which the gas is radially still
supersonic. This layered accretion is described in Tanigawa
et al. (2012) and was also seen by Takasao et al. (2021). In their
simulations, the layer, which is resolved, is much thicker and of
order of 15°. We discuss this in Section 5.2.

With our fiducial resolution of $N_{\theta}=181$, we obtain layered
accretion and outward-directed “backflows” in the layer right
below the shock (Tanigawa et al. 2012; Takasao et al. 2021), in
which $v_{\theta}$ is subsonic but $v_r$ is still supersonic (see Figure 9).
With $N_{\theta}=51$, the gas directly settles vertically to the midplane
instead of performing a “U-turn” in a thin layer beforehand.
However, the fully supersonic part of the flow, especially close
to the planet, is independent of the resolution, and the $H\alpha$
emission also does not depend on the resolution.

Finally, we see that the behavior of the gas at the outer edge
depends somewhat on the resolution. For $N_{\theta}=181$, the
streamlines within about 4° of the midplane flow downward,
while the others are lifted up (Figure 1(b)). At $N_{\theta}=51$, all
streamlines are lifted up. Also, at the outer edge of the CPD, the
streamlines are lifted up in the $N_{\theta}=181$ but not the $N_{\theta}=51$
simulation.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Streamlines for simulations at two $\theta$ resolutions: $\Delta\theta = 0.5^\circ$ (fiducial run; red) and $1.8^\circ$ (black). The density structure (essentially the same for both) is shown in Figure 1. In each panel, the $\theta_{\text{init}} = \theta(r_{\text{max}})$ values of the streamlines are the same for both simulations. Each panel focuses on a different region at a different scale. The flow beneath the CPD shock depends qualitatively on the resolution in the polar direction, but this does not affect our results.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

## Appendix C
### Variations in the Flow Pattern

Ideally, we would be able to wait for a quasi-steady state to establish in the flow at large and small scales and measure the different properties ($\dot{M}_{\text{H}\alpha}$, etc.) from this. In practice, despite months of wall-clock runtime, in the fiducial simulation, a density wave was still traveling out in Phase I (as mentioned in Section 3.2.4). It is associated with the growing CPD outer radius and reflects our setup in which we let the simulation begin without a CPD. This wave somewhat changes the angular distribution of the mass infall close to the freeze radius and thus, in principle, close to the shock radius for Phase II.

We assess how much variation in the H$\alpha$-generating accretion rate could come from this wave. For this, we measure as a function of time the mass flux in the supersonic region within $45^\circ$ of the pole at a distance of $r = 15, 20, 30$, and $50 R_J$ from the planet. Since the flow is smooth, these partial accretion rates will correlate directly with $\dot{M}_{\text{H}\alpha}$, which is not accessible in Phase I because the freeze radius $r_{\text{frz}}$ is farther out than the maximal radius for H$\alpha$ generation. Furthermore, since even after $t = 900 t_{\text{ff,glob}}$, the density wave has not yet reached $r_{\text{max}}$ but rather is still moving out, we look at simulations with different parameters in which the evolution happens more quickly.

We find that in a simulation with identical parameters but a surface density increased by a factor of 10, large-scale oscillations begin around $700 t_{\text{ff,glob}}$. The maximum mass flux is $\dot{M} = 10^{-5} M_J \text{ yr}^{-1}$ at $r = 50 R_J$ and goes smoothly down as a power law to $5 \times 10^{-6} M_J \text{ yr}^{-1}$ at $15 R_J$. The minimum mass flux decreases more steeply from $\dot{M} = 2 \times 10^{-7} M_J \text{ yr}^{-1}$ at $r = 50 R_J$ to $2 \times 10^{-8} M_J \text{ yr}^{-1}$ at $15 R_J$. The total influx at $R_{\text{Hill}}$ is $\dot{M}_{\text{Hill, net}} = 7 \times 10^{-5} M_J \text{ yr}^{-1}$ at all times. Therefore, extrapolating down to a radial distance of $r \approx 2 R_J$, the planet-reaching or H$\alpha$-generating accretion rate is in the range of $0.05 \sim 10^{-5}$ times $\dot{M}_{\text{Hill, net}}$. The minimum value has a considerable uncertainty due to the extrapolation. These are partial rates and are integrated in angle only down to $45^\circ$ from the pole, but the correction down to the CPD height would not

be too large. Assuming that these relative numbers are independent of the surface density and thus also apply to the fiducial run, we would obtain $\dot{M}_{\text{p, direct}}$ or $\dot{M}_{\text{H}\alpha}$ values only up to a factor of $\approx 20$ larger than what we found in the fiducial run (see Figure 4(b)) if we let Phase II begin from a different moment of Phase I. At the other extreme, the partial mass fluxes could be orders of magnitude smaller than what we found.

The upshot of this estimate is that there are transient oscillations, but they will not affect the basic and crucial point that only a fraction, clearly below 100%, of the gas falling onto the CPD can generate emission lines.

## Appendix D
### Temperature and Luminosity Structure

Figure 10 shows the temperature near the planet surface, and Figure 11 shows the velocity, temperature, fluxes, and angular

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Temperature close to the planet. The color scale is capped at 3000 K, but the Zel'dovich spikes—the extremely thin bright regions above the planetary and CPD surface—as well as the innermost hydrostatic parts near $r_{\text{min}}$ reach much higher temperatures. Streamlines for the bolometric radiation are shown (dashed white lines), whereas hydrogen lines, including H$\alpha$, originate from both shocks (see regions with a nongray preshock velocity in Figure 2). The radial segment at $\theta = 77^\circ$ (green) is analyzed in Figure 11 and compared to a vertical segment at $R = 2.5 R_J$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** The 1D cuts along a line of constant angle $\theta = 77^\circ$ (left column; corresponding to the green line segment in Figure 10, just above the CPD surface) and constant cylindrical radius $R = 2.5 R_J$ (right column) in the $N_\theta = 181$ simulation. In the pink regions, the gas is in hydrostatic equilibrium. (a) and (f) Radial $v_r$, total poloidal $v_{pol}$, and freefall velocities $v_{ff,\infty}$. (b) and (g) Gas and radiation temperatures, showing in panel (b) the radial zoning, and with $T_{acc}$. (c) and (h) Reduced flux in $r$ and $\theta$. (d) and (i) Bolometric radiative flux $F_{\text{rad}, r}$ and mechanical energy flux $F_{\text{kin}, r} = 0.5 \rho |v_r|^3$ in the direction $x=r$ or $x=\theta$, corresponding to $T_{\text{acc}}$. (e) and (j) Angular frequency $\Omega = v_\phi/R$ compared to the Keplerian frequency $\Omega_{\text{Kep}}$. Panel (e) also shows $\Omega$ in the midplane (dashed), which is in (rotation-modified) hydrostatic equilibrium in $R$ and $z$ out to the CPD outer edge (not shown). The apparent steps in the right column are a plotting artifact.</figcaption></figure>

frequency along two cuts. In Figure 10, the shocks on the planet and CPD surfaces are clearly visible as Zel'dovich spikes (Zel'dovich & Raizer 1967; see also discussion in Paper II). Thanks to the small cell sizes (see Figure 8), they reach, respectively, $T > 25,000$ and $>4000$ K, off the color scale (capped at 3000 K), but this is resolution-dependent. The true physical peak temperature would be of order $10^4-10^5$ K

(Aoyama et al. 2018). Fortunately, this need not be resolved to follow the radiation transfer correctly (Paper II). The pre- and postshock temperatures, which set more directly the thermal structure of the accretion flow and the settling layers below the shock, are equal and near $T = 1100$ K.

Below the CPD surface shock, the temperature is nearly constant (except close to the midplane), which reflects the low
opacity. Nevertheless, the polar reduced flux,<fnref n="12" />

$$f_{\text{red},\theta} \equiv \frac{|F_{\text{rad},\theta}|}{cE_{\text{rad}}}, \qquad (\text{D1})$$

where $F_{\text{rad}, \theta}$ is the radiation flux in the polar direction, is at most $f_{\text{red}, \theta} \approx 0.1$, while the radial reduced flux $f_{\text{rad}, r}$ goes smoothly from $f_{\text{red}, r} = 0.01-0.1$ near the midplane to $f_{\text{red}, r} = 1$ below, at, and above the CPD surface shock. Thus, the radiation diffuses in the polar direction while also diffusing radially (below the CPD shock) or flowing freely (above it).

The temperature at the shock on the planet surface is $T_{\text{shock}} \approx 1150$ K. This, however, is set mostly by the luminosity below the shock coming from the compression of the gas, namely, the free-streaming “accretion temperature” for an $\eta^{\text{kin}} = 100\%$ shock efficiency (Paper I), given by

$$\sigma_{\text{SB}} T_{\text{acc}}^4 = F_{\text{kin},r} = \frac{1}{2} \rho |v_r|^3, \qquad (\text{D2})$$

where $\sigma_{\text{SB}}$ is the Stefan-Boltzmann constant and is only $T_{\text{acc}} = 615$ K at $\theta = 77^{\circ}$ or $T_{\text{acc}} = 715$ K at the pole. In both cases, this is much smaller than $T_{\text{shock}}$. (This is the limit $\ell = 1$ of Equation (33) in Paper II, while here $\ell \gg 1$, since the downstream luminosity dominates.) In the classical assumption of pure radial infall, the direct infall $\dot{M}_{\text{p, direct}}$ would be predicted to lead to an accretion temperature $T'_{\text{acc}} = (GM_{\text{p}}\dot{M}_{\text{p, direct}}/[4\pi R_{\text{p}}^3\sigma_{\text{SB}}])^{1/4} = 684$ K, ignoring a factor $f_{\text{fill}}/\zeta$ (Zhu 2015), of order unity. As it should, $T'_{\text{acc}}$ lies between the pole and equator values for $T_{\text{acc}}$. However, the pendant to this (from a global simulation point of view) is to implicitly assume that the entire mass flux $\dot{M}_{\text{Hill, net}}$ shocks on the planetary surface, leading to $T_{\text{acc, class}} = 2665$ K, which would dominate the interior luminosity. Neither this radiation temperature $T_{\text{acc, class}}$ nor the corresponding gas temperature in the free-streaming limit $T = T_{\text{acc, class}}/4^{1/4}$ have any relevance in describing the system; the gas falls in more slowly and spreads over a much larger area than assumed by the formula.

On the surface of the CPD at $R = 2.5 R_{\text{J}}$, the temperature is $T = 1070$ K, with the actual $T_{\text{acc}} = 505$ K again much smaller in terms of the radiation fluxes $F \propto T^4$ (the gas and radiation temperatures are equal). Thus, also for the CPD, it is the interior luminosity, not the kinetic energy of the gas, that is responsible for setting the temperature.

In the midplane, there is no shock at the planet surface. Instead, the planet and CPD are connected by a boundary layer (e.g., Hertfelder & Kley 2017; Dong et al. 2021) in which the angular velocity in the midplane peaks somewhat above the Keplerian value $\Omega_{\text{Kep}} = \sqrt{GM_{\text{p}}/r^3}$ before decreasing smoothly to join the boundary condition at $r_{\text{min}}$ (Figure 11(e)). This region will not be studied further here. At least at $R = 2.5 R_{\text{J}}$, the whole vertical extent of the CPD is in Keplerian rotation; in the pink region in Figure 11(j), $\Omega = \Omega_{\text{Kep}}$. The boundary layer leads to a higher temperature close to the midplane, but only slightly so.

Away from the CPD (for $r \gtrsim 100 R_{\text{J}}$), the temperature distribution is independent of polar angle, with temperatures below 100 K. In that regime, the dust opacity $\kappa_{\text{dust}} \approx 0.01 \text{ cm}^2 \text{ g}^{-1}$ dominates by 3–4 dex over the gas opacity,

even for our choice of a low $f_{\text{d/g}} = 10^{-4}$. The radial Rosseland optical depth from $r_{\text{max}}$ to the shock is $\Delta\tau_{\text{Ross}} \sim 3 \times 10^{-3}$ along the pole, roughly a factor of 2 higher on a path just above the CPD, and $\Delta\tau_{\text{Ross}} = 10^{-3}$ in the midplane down to the outer edge of the CPD. The low overall optical depth reflects the low $f_{\text{d/g}}$ and modest mass inflow into the Hill sphere $\dot{M}_{\text{Hill, net}} \sim 10^{-5} M_{\text{J}} \text{ yr}^{-1}$ (Table 2).

In this particular example, most of the bolometric flux reaching the observer is coming from the interior of the planet and from the CPD itself. These fluxes do not come from the immediate conversion of kinetic energy but rather from the cooling of the hydrostatic regions below the shocks. This, in turn, depends on the accretion history. In a given simulation, this history is set by the numerical approach (here the two-phase system we used, which spans several freefall timescales) and in general by the variation of the accretion rate over formation timescales of order 1 Myr.

In Section 5.1, we compared the $H\alpha$ flux that we predict to the observed one for PDS 70 b. Assuming a roughly linear scaling of the $H\alpha$ luminosity with the mass inflow rate, the latter would need to be 7–100 times larger than in our simulation in order to match the observed $L_{\text{H}\alpha}$. From Equation (D2), this would imply $T_{\text{acc}} \approx 1000-1900$ K at the planet’s surface near the CPD and up to $T_{\text{acc}} \approx 2300$ K at the pole. In this case, the accretion luminosity from the shock would likely dominate the temperature structure, and the highest accretion rates would be more challenging to reconcile with the constraints from Wang et al. (2021) on $T_{\text{eff}}$ from the K-band spectral shape. Details such as the viewing geometry or complex radiation transfer effects could, however, play an important role. Next-generation spectroscopic observations would be helpful to develop a robust and self-consistent picture.

## ORCID iDs

Gabriel-Dominique Marleau https://orcid.org/0000-0002-2919-7500

Rolf Kuiper https://orcid.org/0000-0003-2309-8963

Christoph Mordasini https://orcid.org/0000-0002-1013-2811

## References

Adams, F. C., & Batygin, K. 2022, ApJ, 934, 111

Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84

Aoyama, Y., Marleau, G.-D., Ikoma, M., & Mordasini, C. 2021, ApJL, 917, L30

Aoyama, Y., Marleau, G.-D., Mordasini, C., & Ikoma, M. 2020, arXiv:2011.06608

Asensio-Torres, R., Henning, T., Cantalloube, F., et al. 2021, A&A, 652, A101

Ayliffe, B. A., & Bate, M. R. 2009a, MNRAS, 393, 49

Ayliffe, B. A., & Bate, M. R. 2009b, MNRAS, 397, 657

Ayliffe, B. A., & Bate, M. R. 2012, MNRAS, 427, 2597

Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJL, 884, L41

Bailey, A., Stone, J. M., & Fung, J. 2021, ApJ, 915, 113

Batygin, K. 2018, AJ, 155, 178

Berardo, D., Cumming, A., & Marleau, G.-D. 2017, ApJ, 834, 149

Béthune, W., & Rafikov, R. R. 2019a, MNRAS, 487, 2319

Béthune, W., & Rafikov, R. R. 2019b, MNRAS, 488, 2365

Betti, S. K., Follette, K. B., Ward-Duong, K., et al. 2022a, ApJL, 935, L18

Betti, S. K., Follette, K. B., Ward-Duong, K., et al. 2022b, ApJL, 941, L20

Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icar, 143, 2

Bonse, M. J., Garvin, E. O., Gebhard, T. D., et al. 2023, arXiv:2303.12030

Brittain, S. D., Najita, J. R., Dong, R., & Zhu, Z. 2020, ApJ, 895, 48

Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, NatAs, 2, 138

Bryan, M. L., Ginzburg, S., Chiang, E., et al. 2020, ApJ, 905, 37

Chachan, Y., Lee, E. J., & Knutson, H. A. 2021, ApJ, 919, 63

<footnote n="12">¹² The reduced flux, or “streaming factor” (Kley 1989), measures the extent to which radiation is diffusing ($f_{\text{red}} \to 0$) or freely streaming ($f_{\text{red}} \to 1$).</footnote>
Chazelas, B., Lovis, C., Blind, N., et al. 2020, Proc. SPIE, 11448, 1144875

Chen, Z., & Bai, X. 2022, ApJL, 925, L14

Christensen, U. R., Holzwarth, V., & Reiners, A. 2009, Natur, 457, 167

Choksi, N., Chiang, E., Fung, J., & Zhu, Z. 2023, arXiv:2305.01684

Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662

Close, L. M. 2020, AJ, 160, 221

Cugno, G., Quanz, S. P., Hunziker, S., et al. 2019, A&A, 622, A156

Currie, T., Lawson, K., Schneider, G., et al. 2022, NatAs, 6, 751

Delorme, J.-R., Jovanovic, N., Echeverri, D., et al. 2021, JATIS, 7, 035006

Dong, J., Jiang, Y.-F., & Armitage, P. J. 2021, ApJ, 921, 54

Dražkowska, J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H. 2019, ApJ, 885, 91

Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A69

Eriksson, S. C., Asensio Torres, R., Janson, M., et al. 2020, A&A, 638, L6

Follette, K. B., Close, L. M., Males, J. R., et al. 2023, AJ, 165, 225

Fu, Z., Huang, S., & Yu, C. 2023, ApJ, 945, 165

Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152

Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125

Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, NatAs, 3, 749

Haffert, S. Y., Males, J. R., Close, L., et al. 2021, Proc. SPIE, 11823, 1182306

Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, ApJ, 475, 770

Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135

Hasegawa, Y., Kanagawa, K. D., & Turner, N. J. 2021, ApJ, 923, 27

Hertfelder, M., & Kley, W. 2017, A&A, 605, A24

Hill, G. W. 1878, AmJM, 1, 5

Huélamo, N., Chauvin, G., Mendigutía, I., et al. 2022, A&A, 668, A138

Isella, A., & Natta, A. 2005, A&A, 438, 899

Kanagawa, K. D., Tanaka, H., Muto, T., & Tanigawa, T. 2017, PASJ, 69, 97

Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140

Karlin, S. M., Panić, O., & van Loo, S. 2023, MNRAS, 520, 1258

Katarzyński, K., Gawroński, M., & Goździewski, K. 2016, MNRAS, 461, 929

Keith, S. L., & Wardle, M. 2014, MNRAS, 440, 89

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Kley, W. 1989, A&A, 208, 98

Kley, W. 1998, A&A, 338, L37

Korycansky, D. G., & Papaloizou, J. C. B. 1996, ApJS, 105, 181

Krapp, L., Kratter, K. M., & Youdin, A. N. 2022, ApJ, 928, 156

Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010, A&A, 511, A81

Kuiper, R., Yorke, H. W., & Mignone, A. 2020, ApJS, 250, 13

Kurokawa, H., & Tanigawa, T. 2018, MNRAS, 479, 635

Lambrechts, M., & Lega, E. 2017, A&A, 606, A146

Lambrechts, M., Lega, E., Nelson, R. P., Crida, A., & Morbidelli, A. 2019, A&A, 630, A82

Lovelace, R. V. E., Covey, K. R., & Lloyd, J. P. 2011, AJ, 141, 51

Machida, M. N., Kokubo, E., Inutsuka, S.-i., & Matsumoto, T. 2008, ApJ, 685, 1220

Maeda, N., Ohtsuki, K., Tanigawa, T., Machida, M. N., & Suetsugu, R. 2022, ApJ, 935, 56

Mai, C., Desch, S. J., Kuiper, R., Marleau, G.-D., & Dullemond, C. 2020, ApJ, 899, 54

Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, A91

Marleau, G.-D., Aoyama, Y., Kuiper, R., et al. 2022, A&A, 657, A38

Marleau, G.-D., Klahr, H., Kuiper, R., & Mordasini, C. 2017, ApJ, 836, 221, (Paper I)

Marleau, G.-D., Mordasini, C., & Kuiper, R. 2019, ApJ, 881, 144, (Paper II)

Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541

Marois, C., Lafrenière, D., Macintosh, B., & Doyon, R. 2008, ApJ, 673, 647

Mendoza, S., Tejeda, E., & Nagel, E. 2009, MNRAS, 393, 579

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7

Moldenhauer, T. W., Kuiper, R., Kley, W., & Ormel, C. W. 2021, A&A, 646, L11

Moldenhauer, T. W., Kuiper, R., Kley, W., & Ormel, C. W. 2022, A&A, 661, A142

Mordasini, C. 2013, A&A, 558, A113

Mordasini, C., Alibert, Y., Georgy, C., et al. 2012a, A&A, 547, A112

Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012b, A&A, 547, A111

Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2

Mungan, C. E. 2009, PhTea, 47, 502

Nelson, R. P., Lega, E., & Morbidelli, A. 2023, A&A, 670, A113

Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13

Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10

Pringle, J. E. 1981, ARA&A, 19, 137

Ringqvist, S. C., Viswanath, G., Aoyama, Y., et al. 2023, A&A, 669, L12

Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2002, ApJ, 578, 420

Sanchis, E., Picogna, G., Ercolano, B., Testi, L., & Rosotti, G. 2020, MNRAS, 492, 3440

Sanghi, A., Zhou, Y., & Bowler, B. P. 2022, AJ, 163, 119

Schulik, M., Johansen, A., Bitsch, B., & Lega, E. 2019, A&A, 632, A118

Schulik, M., Johansen, A., Bitsch, B., Lega, E., & Lambrechts, M. 2020, A&A, 642, A187

Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Stolker, T., Marleau, G. D., Cugno, G., et al. 2020, A&A, 644, A13

Szulágyi, J. 2017, ApJ, 842, 103

Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1

Szulágyi, J., Dullemond, C. P., Pohl, A., & Quanz, S. P. 2019, MNRAS, 487, 1248

Szulágyi, J., & Ercolano, B. 2020, ApJ, 902, 126

Takasao, S., Aoyama, Y., & Ikoma, M. 2021, ApJ, 921, 10

Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47

Thanathibodee, T., Calvet, N., Bae, J., Muzerolle, J., & Hernández, R. F. 2019, ApJ, 885, 94

Thommes, E. W., Matsumura, S., & Rasio, F. A. 2008, Sci, 321, 814

Toci, C., Lodato, G., Christiaens, V., et al. 2020, MNRAS, 499, 2015

Ulrich, R. K. 1976, ApJ, 210, 377

Uyama, T., Xie, C., Aoyama, Y., et al. 2021, AJ, 162, 214

Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72

Wagner, K., Apai, D., & Kratter, K. M. 2019, ApJ, 877, 46

Wagner, K., Follette, K. B., Close, L. M., et al. 2018, ApJL, 863, L8

Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148

Xie, C., Haffert, S. Y., de Boer, J., et al. 2020, A&A, 644, A149

Zel'dovich, Y. B., & Raizer, Y. P. 1967, Physics of Shock Waves and High-temperature Hydrodynamic Phenomena (New York: Academic Press)

Žerjal, M., Ireland, M. J., Crundall, T. D., Krumholz, M. R., & Rains, A. D. 2023, MNRAS, 519, 3992

Zhou, Y., Bowler, B. P., Wagner, K. R., et al. 2021, AJ, 161, 244

Zhou, Y., Sanghi, A., Bowler, B. P., et al. 2022, ApJL, 934, L13

Zhu, W., & Dong, S. 2021, ARA&A, 59, 291

Zhu, Z. 2015, ApJ, 799, 16

Zurlo, A., Cugno, G., Montesinos, M., et al. 2020, A&A, 633, A119