# Spectral appearance of the planetary-surface accretion shock:
Global spectra and hydrogen-line profiles and luminosities

Yuhiko Aoyama Department of Earth and Planetary Science, University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Institute for Advanced Study,
Tsinghua University, Beijing 100084, People’s Republic of China  Department of
Astronomy, Tsinghua University, Beijing 100084, People’s Republic of China
Gabriel-Dominique Marleau Institut für Astronomie und Astrophysik, Universität
Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany  Physikalisches
Institut, Universität Bern, Gesellschaftsstr. 6, 3012 Bern, Switzerland  Max-
Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
Christoph Mordasini Physikalisches Institut, Universität Bern,
Gesellschaftsstr. 6, 3012 Bern, Switzerland  Masahiro Ikoma Department of
Earth and Planetary Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan  Division of Science, National Astronomical Observatory
of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

(Received —; Revised —; Accepted —)

###### Abstract

Hydrogen-line emission from an accretion shock has recently been observed at
planetary-mass objects. Our previous work predicted the shock spectrum and
luminosity for a shock on the circumplanetary disc. We extend this to the
planet-surface shock. We calculate the global spectral energy distribution
(SED) of accreting planets by combining our model emission spectra with
photospheric SEDs, and predict the line-integrated flux for several hydrogen
lines, especially H$\alpha$, but also H$\beta$, Pa$\alpha$, Pa$\beta$,
Pa$\gamma$, Br$\alpha$, and Br$\gamma$. We apply our non-equilibrium emission
model to the surface accretion shock for a wide range of accretion rates
$\dot{M}$ and masses $M_{\mathrm{p}}$. Fits to formation calculations provide
radii and effective temperatures. Extinction by the surrounding material is
neglected, which is arguably often relevant. We find that the line luminosity
increases monotonically with $\dot{M}$ and $M_{\mathrm{p}}$, depending mostly
on $\dot{M}$ and weakly on $M_{\mathrm{p}}$ for the relevant range of
parameters. The Lyman, Balmer, and Paschen continua can exceed the
photosphere. The H$\beta$ line is fainter by 0–1 dex than H$\alpha$, whereas
other lines are weaker (by $\sim 1$–3 dex). Shocks on the planet or the CPD
surface are distinguishable at very high spectral resolution, but the planet
surface shock likely dominates if both are present. Applied to recent non-
detections of H$\alpha$, our models imply looser constraints on the $\dot{M}$
of putative planets than when extrapolating fits from the stellar regime.
These hydrogen-line luminosity predictions are useful for interpreting
(non-)detections of accreting planets.

Exoplanet formation (492); Accretion (14); Shocks (2086); Hydrogen lines [H i
line emission] (690); H alpha photometry (691); Direct imaging (387); High
resolution spectroscopy (2096)

††journal: ApJ

##  1 Introduction

Recent instrumental improvement have enabled the observation of forming
planets (e.g. Kraus & Ireland, 2012; Quanz et al., 2013; Currie et al., 2015;
Wagner et al., 2018). Because it is expected to yield information on how
planets grow, the detection of H $\alpha$ is particularly important (Sallum et
al., 2015; Wagner et al., 2018; Haffert et al., 2019; Cugno et al., 2019;
Zurlo et al., 2020; Xie et al., 2020).

In the context of forming low-mass protostars (Classical T Tauri stars: CTTS),
H $\alpha$ is known as an indicator of accretion and used to estimate mass
accretion rate (e.g. Gullbring et al., 1998). The H $\alpha$ from CTTS is
brighter than the photospheric continuum by a few mag and has large width
($\gtrsim 200~{}\mathrm{km}\,\mathrm{s^{-1}}$). The magnetospheric accretion
model (Uchida & Shibata, 1984; Königl, 1991) can explain these characteristic
features when the accretion funnel is hot enough to emit H $\alpha$ (e.g.
Hartmann et al., 1994; Muzerolle et al., 2001). Furthermore, the line-
integrated luminosity ($L_{\mathrm{H}\,\alpha}$) or the spectral width
(${\Delta\lambda_{\textnormal{H}\,\alpha}}$) of the H $\alpha$ line shows a
correlation to the mass-accretion rate (or accretion luminosity,
$L_{\mathrm{acc}}$) estimated with modeling of continuum emission (e.g.
Valenti et al., 1993; Calvet & Gullbring, 1998). Therefore, H $\alpha$ is used
to estimate the accretion rate for protostars that are too far for their
continuum emission to be observable (e.g. Gullbring et al., 1998; Herczeg &
Hillenbrand, 2008; Fang et al., 2009; Rigliaco et al., 2012; Alcalá et al.,
2014, 2017; Natta et al., 2004). Similar statements hold for further hydrogen
lines from the Balmer, Paschen, Brackett, or other series.

As for the stellar case, an H $\alpha$ excess was reported from (candidate)
protoplanets, and the observed luminosity was used to estimate the accretion
rate by applying the results of the CTTS observations (Sallum et al., 2015;
Wagner et al., 2018; Haffert et al., 2019). However, there is no guarantee
that relationships between $L_{\mathrm{acc}}$ and $L_{\mathrm{H}\,\alpha}$ or
${\Delta\lambda_{\textnormal{H}\,\alpha}}$ given in CTTS are valid for
protoplanets. Thanathibodee et al. (2019) applied the stellar H $\alpha$
emission model of Muzerolle et al. (2001) to a planetary-mass object (PDS 70 b
(catalog )) and argued the $L_{\mathrm{H}\,\alpha}$–$L_{\mathrm{acc}}$
relationship shows a different trend from that of protostars (Ingleby et al.,
2013; Rigliaco et al., 2012; see also Szulágyi & Ercolano, 2020 and the
discussion of their work in Aoyama et al. submitted).

Planetary gas accretion is qualitatively different from the stellar one in
some points. An important characteristic feature is that protoplanets and
their surrounding gaseous disk (circum-planetary disk, CPD) are embedded in
the stellar surrounding disk (protoplanetary disk, PPD). On the way of gas
accreting towards the protoplanet, the gas preferentially enters the planetary
gravitational sphere in high altitudes above the disk midplane (e.g. Tanigawa
et al., 2012). When the gas falling from the PPD to CPD vertically hits the
CPD surface, it yields a strong shock, which can be hot enough to emit H
$\alpha$ (Szulágyi & Mordasini, 2017). Aoyama et al. (2018) constructed a
model of shock-heated gas with cooling, chemical reactions, and radiative
transfer, estimated hydrogen line luminosity depending on the gas velocity and
density, and estimated the $L_{\mathrm{H}\,\alpha}$ depending on the shock
properties.

On the other hand, the magnetospheric accretion may occur even in the
planetary accretion, bringing about a strong shock also on the planetary
surface. If protoplanets have dipole magnetic fields strong enough to control
the gas dynamics, vertical accretion can occur directly onto the planetary
surface (Batygin, 2018). While the accretion shock on the CTTS surface is too
strong and makes gas too hot to emit H $\alpha$ (see e.g., Hartmann et al.,
2016), the weak gravity of protoplanets leads to moderate free-fall velocity
($\sim 100~{}\mathrm{km}\,\mathrm{s^{-1}}$) and to emitting a significant
amount of H $\alpha$. In contrast, in the CPD surface shock model, only a
small fraction ($\lesssim 1\%$) can contribute to the H $\alpha$ emission,
because most gas hits the CPD far from the planet (Aoyama et al., 2018). Also,
in the magnetospheric accretion-funnel model, the heating mechanism is still
an open question (Muzerolle et al., 2001). Therefore, the gas in the accretion
funnel could be too cool to emit H $\alpha$, perhaps especially for
protoplanets not much more massive than Jupiter.

This motivates us, in this study, to model the hydrogen line emission coming
from the planetary surface shock, considering a wide range of parameters. We
focus on H $\alpha$ first and then explore other hydrogen line emission. We
combine these results with models of the photospheric emission and discuss
when the shock lines are visible above the photosphere emission. Note that
part of the planetary surface shock model presented here was used already in
Aoyama & Ikoma (2019) for the case of PDS 70 b and c. A more extensive
investigation is done in this study.

The paper is organized as follows: In Section 2 we discuss the properties of
the planetary-surface shock and of the planets and review our numerical shock
model, which was introduced in Aoyama et al. (2018). In Section 3 we present
emission spectra of accreting gas giants for a large grid of models, before
applying this in Section 4 to a few objects, especially to their detection at
H $\alpha$. In Section 5 we explore other observational aspects, including
line strengths for lines other than H $\alpha$ and the possibility of breaking
some degeneracies. Finally, we present a critical discussion in Section 6
before summarizing in Section 7. The appendices present further material: a
discussion of our approach compared to Storey & Hummer (1995) (Appendix B),
the inverse relationship between the shock-microphysical and planet-formation
parameters (Appendix C), a map of the H $\alpha$ luminosity for the cold-start
radius fits (Appendix D), and the calculation of the H $\alpha$ luminosity in
Wagner et al. (2018) (Appendix E).

##  2 Description of the combined model

We model the spectral energy distribution (SED) of an accreting gas giant with
a surface accretion shock. The radiation from the accreting gas giant is
mainly composed of two components, namely the photospheric radiation (Section
2.2) and the shock excess (Section 2.3). We assume that the two components can
be computed separately, i.e., that the layers heated by the shock do not
affect significantly the rest of the emission.

The input parameters for our combined spectra of the accretion shock and the
photosphere are the following five: mass accretion rate $\dot{M}$, planet mass
$M_{\mathrm{p}}$, planet radius $R_{\mathrm{p}}$, filling factor
$f_{\mathrm{fill}}$ of the shock on the planet surface, and photospheric
effective temperature ${T_{\textnormal{eff}}}$. However, taking them as
independent would result in an impractically large parameter space and may
lead to unlikely combinations (e.g., small radius and mass but large
luminosity). Therefore, in this study, only ($\dot{M},M_{\mathrm{p}}$) will be
considered as free parameters via the modeling described in Section 2.1. We
consider here the case that the planetary emission (photosphere and shock) is
not extincted, and detail in Section 2.4 when this is relevant.

###  2.1 Fitting of planetary properties

For convenience, the planet radius and effective temperature are derived from
a specific detailed planet formation and structure model. The planet radius
$R_{\mathrm{p}}$ is fitted as a function of ($\dot{M},\,M_{\mathrm{p}}$) by
using the Bern model (Alibert et al., 2005; Mordasini et al., 2012b, a, 2015,
2017; Emsenhuber et al., 2020a, b; Schlecker et al., 2020). We use a constant
intrinsic temperature of $1000$ K and the accretion heating to predict the
effective temperature ${T_{\textnormal{eff}}}$. For consistency, we assume, in
the incoming mechanical energy at the shock, that the residual that is not
radiated from the shock emission heats the photosphere. More details are given
in Section A.1 (see Equation A5).

###  2.2 Photospheric emission

For the photospheric radiation model, we use the CIFIST2011_2015 BT-Settl
models, which calculated spherical radiative transfer in atmospheres with
solar metallicity111From https://phoenix.ens-lyon.fr/Grids/BT-
Settl/CIFIST2011_2015/. (Allard et al., 2012; Baraffe et al., 2015). This
requires the effective temperature, surface gravity, and the emitting area.
They are derived from
(${T_{\textnormal{eff}}},\,M_{\mathrm{p}},\,R_{\mathrm{p}},\,f_{\mathrm{fill}}$).

The BT-Settl model simulates the photospheric emission from isolated objects.
The accretion heats the top of the atmosphere and, in general, will change the
temperature structure. Since the detailed absorption feature could highly
depend on the temperature structure in the upper layer, they are less reliable
for accreting objects. However, this model can show how bright the shock
emission is relative to the photospheric emission, i.e., the detectability of
shock emission.

In this study, we focus on the shock-heated gas on the planetary surface. We
treat only the emission from the photosphere and shock-heated gas but not from
the CPD, whose temperature is lower than those of the photosphere and the
shock. Continuum emission from a (simplified) CPD model has been calculated in
Zhu (2015), Eisner (2015), and Szulágyi et al. (2019), and the line emission
from the shock on the CPD has been calculated in Aoyama et al. (2018).

###  2.3 Shock emission

The shock excess is calculated from the 1D radiation-hydrodynamic model
developed by Aoyama et al. (2018) and Aoyama & Ikoma (2019), which is outlined
in Section A.2. This model predicts the shock emission flux from the gas
velocity $v_{0}$ and number density $n_{0}$ just before the shock. These shock
parameters and the emitting area can be derived with
($\dot{M},\,M_{\mathrm{p}},\,R_{\mathrm{p}},\,f_{\mathrm{fill}}$), assuming
$v_{0}$ to be the free-fall velocity. The details are described in Section
A.2. For simplicity, we use $f_{\mathrm{fill}}=1$ because $f_{\mathrm{fill}}$
hardly changes the line flux (see also Section 6.2). A discussion of how the
accretion geometry sets $f_{\mathrm{fill}}$ is given in Marleau et al. (subm.,
see their Figure 1).

In this study, we treat only the shock on the planetary surface rather than
that on CPD surface, which Aoyama et al. (2018) studied. We compare the
emission from both in Section 5.2.

Also, we contrast our model of the shock emission to Storey & Hummer (1995) in
Appendix B. In short, for the shock heated gas, collisional excitation from
the ground state is a crucial element that we include, while SH95 did not.

###  2.4 Neglecting extinction

In this work, extinction by material between the shock surface (the planet)
and the observer is not considered. Several components contribute to it. In
principle, the contribution from the interstellar medium (ISM) can be
determined for a given source, from the stellar spectrum or by statistical
tools such as Stilism of Lallement et al. (2019). Thus this component is
relatively easily accounted for. To what extent the gas or dust surrounding a
forming planet may weaken the shock signal is an important question, as the
recent observational results and theoretical modeling in Hashimoto et al.
(2020), Stolker et al. (2020a), and Sanchis et al. (2020) highlight. However,
considering extinction adds an entire level of complexity and brings many
uncertainties, in particular concerning the radiative transfer geometry and
the dust opacity. Therefore, we deal with extinction by the gas and the dust
in a dedicated paper (Marleau et al., subm.).

Nevertheless, the extinction-free case is relevant in itself. As we show in
Marleau et al. (subm.), an accretion flow free of extinction at H $\alpha$ is
a plausible assumption for a wide range of accretion rates and masses. We look
at this in detail but, heuristically, the transition disk gaps in which
planets are found are usually dust-free (Close, 2020), and gas cooler than a
few thousand kelvin can be optically thin for non-Lyman-series hydrogen lines.
Also, while Hashimoto et al. (2020) derived for PDS 70 b an extinction of
$>2.0$ mag at H $\alpha$, we should recall that this is for a single object
(PDS 70 b), and that this estimate depends on the wide spectral width of the
observed H $\alpha$ line (see Section 4.1), which can be overestimated due to
the finite instrumental resolution (Thanathibodee et al., 2019). More
generally, it is conceivable that for some accretion and viewing geometries
the H $\alpha$ produced at the shock could leave the system without passing
through any absorbing material that could be present. For these reasons, it
seems sensible to treat the extincted case separately.

##  3 Theoretical spectra of forming gas giants

We now turn to results from the methods described above. We look first at one
representative example in detail (Section 3.1; Aoyama et al. 2018 showed three
other cases for the CPD case) and then survey a large part of the relevant
$(\dot{M},M_{\mathrm{p}})$ parameter space (Section 3.2).

###  3.1 One example

####  3.1.1 Postshock structure

Figure 1: Postshock flow structure, beginning immediately after the
hydrodynamical jump, for $v_{0}=100~{}\mathrm{km}\,\mathrm{s^{-1}}$ and
$n_{0}=10^{11}$ cm-3. This corresponds for example to
$\dot{M}=10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$,
$M_{\mathrm{p}}=5~{}M_{\mathrm{J}}$, and $R_{\mathrm{p}}=1.7~{}R_{\mathrm{J}}$
with $f_{\mathrm{fill}}=1$. The left (right) axes are the depth below (time
elapsed after) the shock surface $\Delta z$ ($\Delta t$). (a): Temperature $T$
(red line, bottom axis) and pressure $P$ (blue line, top linear axis). The
preshock $T$ (not shown) is $T_{0}\approx 1190$ K. The $P$ profile changes
inversely to $T$ because the density increases faster than $T$ drops. (b):
Cooling rates of H collisional excitation and ionization (red line),
dissociation of H2 (black), and OH (brown) and $\mathrm{H_{2}O}$ (blue)
rotational line emission. Throughout the simulation, molecules hardly affect
the cooling because they are minor relative to neutral hydrogen. (c): Number
density relative to $n_{0}$ (see Section A.2.1) for $\mathrm{H}_{2}$ (purple),
$\mathrm{H}^{n=1}$ (black), $\mathrm{H}^{n=2}$ (red), $\mathrm{H}^{n=3}$
(blue), $\mathrm{H}^{+}$ (orange), and $\mathrm{e}^{-}$ (green), with $n$ the
principal quantum number. Thin black lines show $n\geqslant 4$ states. At
intermediate depths, the H ionization fraction approaches unity. (d): Upward
energy flux of Ly $\alpha$ (black line), Ba $\alpha$ (H $\alpha$; red), Pa
$\alpha$ (blue), Pa $\beta$ (green), Br $\gamma$ (brown), recombination
continua (thin dotted; Lyman (black), Balmer (red), Paschen (blue), Blackett
(orange), Pfund (green), Humphreys (pink), and others (purple)), and other
lines (thin black)). At intermediate depths, the line fluxes hardly change due
to the optical thinness of the gas. The deep region in which the fluxes
originate is well resolved by steps of 10 % in $T$, but the details are not
seen easily due to the logarithmic scale.

The postshock structure and hydrogen line emission were detailed by Aoyama et
al. (2018). Although the results shown here are basically the same as theirs,
we review their findings in this subsection for the reader’s convenience.
Also, the input parameters are chosen to be appropriate for the detected
planet PDS 70 b. The $v_{0}$ and $n_{0}$ are higher than in Aoyama et al.
(2018), where we focused on the CPD surface shock rather than the planetary
surface.

In Figure 1 we show one example of the postshock structure for
$v_{0}=100~{}\mathrm{km}\,\mathrm{s^{-1}}$ and $n_{0}=10^{11}$ cm-3. This
corresponds to, for example,
$\dot{M}=10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$,
$M_{\mathrm{p}}=5~{}M_{\mathrm{J}}$, and $R_{\mathrm{p}}=1.7~{}R_{\mathrm{J}}$
with $f_{\mathrm{fill}}=1$, assuming the radius fit to the warm population
(see Figure 12 and Section A.1.1). Following Equation (A19), the preshock
temperature is $T_{0}\approx 1200$ K, and with $\rho_{0}=2.3\times
10^{-13}~{}$g cm-3 this implies $\mu=2.3$ and $\mathcal{M}\approx 45$. The
hydrodynamic shock heats the gas to $T_{1}\approx 4\times 10^{5}$ K (panel
(a)). The gas density (not shown) increases as the temperature drops. Although
the gas cools, the gas pressure increases slightly, by only 30 % at the end of
this simulation, because of the density enhancement of compression. Notice
that the pressure gets low once around the depth of $2\times 10^{3}$ cm
($10^{-3}$ s) because H2 dissociation results in expansion.

Immediately after the shock, the dissociation of H2 is the main process
responsible for the cooling of the gas. However, it can bring down the
temperature only by a small amount before the collisional excitation of the
atomic hydrogen and its ionization take over some $\Delta t\sim 10^{-3}$ s
after the shock (Figure 1b). The excitation and dissociation dominate until
the end of the simulation when $T$ reaches $10^{4}$ K. Throughout the
simulation, molecules hardly affect the cooling because they are minor
relative to neutral hydrogen.

Figure 1c shows the abundance of each form of hydrogen relative to all
hydrogen protons. At the preshock temperature $T_{0}\approx 1200$ K and thus
at the beginning of the evolution in the postshock region, the molecular form
H2 (purple curve) is the most abundant, orders of magnitude more so than the
ground state of atomic hydrogen Hn=1 (black) or ionized hydrogen H+ (dark
yellow). This is why immediately after the shock hydrogen dissociation is the
most important process. Near $\Delta t\approx 10^{-3}$ s, neutral hydrogen
dominates, but very quickly, by $\Delta t\approx 2\times 10^{-3}$ s ($\Delta
z\approx 5000$ cm), the ionization fraction has nearly reached unity.
Excitation from the ground state to the first excited state $n=2$ nevertheless
proceeds simultaneously. The main processes increasing and decreasing the
$n=3$ population is collisional excitation from $n\leqslant 2$ to $n=3$ and
radiative de-excitation from $n=3$ to $n<3$. The $n=2$ level is special
because Ly $\alpha$ ($n=2\rightarrow 1$) is optically thick at this location
in this example, which prevents radiative de-excitation. If the gas were even
denser, also H $\alpha$ would become optically thick, which would prevent
radiative de-excitation and lead to a higher $n=3$ population. At depths
$\Delta z\gtrsim 1\times 10^{7}$ cm, the dropping temperature and the longer
timescale let the hydrogen recombine and the electrons fall back down to lower
levels.

As Figure 1d reveals, the (potentially observable) hydrogen lines originate
mainly in the deepest part, at $\Delta z\approx 2\times 10^{7}$ cm (or $\Delta
t\approx 30$ s), from the de-excitation of the electrons. This region seems
narrow on a logarithmic scale but we recall that the adaptive step size in
time or space is set by a temperature criterion that ensures sufficient
resolution (Aoyama et al., 2018).

Above this region, i.e., for most of the postshock flow as visible in Figure
1, the line fluxes remain approximately constant, with only small modulations
that can be related to the cooling processes shown in panel (b). However, the
shallower region plays an important role in the “wing” of the spectral profile
because of its hot temperature. On the other hand, while the emission from the
deeper region carries most of the energy, its Doppler width is only a few tens
of $\mathrm{km}\,\mathrm{s^{-1}}$ corresponding the temperature of a few
$10^{4}$ K, with hardly a dependence on $v_{0}$.

In general, to zeroth order, when the gas is optically thin and all the energy
is in hydrogen lines, a half goes outward as Ly $\alpha$ and the other half
goes inward also as Ly $\alpha$. In this example, at the shock surface, the
upward-travelling Lyman-$\alpha$ flux represents around 76 % of the incoming
(mostly kinetic) energy, and H $\alpha$ carries only around 1 %. The other
part of the energy influx travels downward, towards the photosphere (see also
Section A.1.2 for more precise fractions).

In our models, we currently do not include cooling from He or metal lines but
this ultimately does not matter. In most regions, hydrogen lines are almost
the only coolant, so that when the abundance of neutral hydrogen becomes low
enough, cooling by hydrogen becomes inefficient. In Figure 1, this is between
$\Delta z\approx 8\times 10^{3}$ and $5\times 10^{7}$ cm. This leads to a
plateau in the temperature, which ends where the hydrogen recombines. In that
temperature region (at $T\sim 10^{5}$ K), cooling by lines of ionic C, O, and
He (specifically, the He i $\lambda 16404$ line) or other metals lines would
be more important (see Figure 3 in Gnat & Ferland, 2012) so that there would
not be a temperature plateau. Indeed, while the ionisation of C and O is
included in the chemistry subroutine, it is not included in the radiation
transfer, and neither is the cooling by lines of C and O in the energy
equation. For helium, the ionisation is included in the energy equation but,
also here, the lines are not. Also, note that for the case presented in Figure
2 of Aoyama et al. (2018), with
$(v_{0}=40~{}\mathrm{km}\,\mathrm{s^{-1}},n_{0}=10^{11}~{}\mathrm{cm}^{-3})$,
in the early parts of the flow the electron abundance is higher than the H+
abundance. These electrons are coming from ionized helium.

However, the gas in the region of the temperature plateau only contributes to
the recombination continuum but not to the hydrogen lines. Therefore, even
including helium or metal lines (and thus changing the temperature structure
of that region) would not modify the strength of the hydrogen lines.

####  3.1.2 Radiative properties

Figure 2: Spectral energy distribution of an accreting gas giant, with
effective temperature ${T_{\textnormal{eff}}}=1400$ K and surface gravity
$\log g=4.0$ (cm s-2). This corresponds for example to a
$M_{\mathrm{p}}=8~{}M_{\mathrm{J}}$ planet accreting at
$\dot{M}=10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ with
$f_{\mathrm{fill}}=1$ in the hot-start population. The flux is for a distance
of 150 pc to the system. The black line shows the pure photospheric radiation
(BT-Settl; Allard et al., 2012) and the red shows the photospheric radiation
with the shock excess. The pale gray line is a blackbody at
${T_{\textnormal{eff}}}=1400$ K. The continua are plotted darker to indicate
that they are less certain. Only the Lyman, Balmer, and Paschen continua are
brighter than the photosphere. Extinction via material between the planet
surface and the observer is not taken into account. UV photons will be
strongly affected, as the gray hatch indicates.

Figure 2 shows the entire corresponding SED, including the contribution from
the photosphere. The ${T_{\textnormal{eff}}}$ is determined by the accretion
heating with a constant intrinsic temperature of 1000 K (see Section A.1),
leading to ${T_{\textnormal{eff}}}=1400$ K. Wherever the photospheric emission
is brighter than the blackbody (gray curve), for instance around
$\lambda\approx 10^{4}$ Å, the photons are coming from layers higher than the
photosphere. In those regions, the actual temperature structure could be
different from the one in BT-Settl due to the shock heating, but this is
difficult to estimate without detailed modeling. For optical or longer
wavelengths ($\lambda\gtrsim 4000$ Å), the radiation is dominated by the
photospheric contribution (black line) except for some hydrogen lines (red
peaks).

On the other hand, at shorter wavelength, the dominant component is Ly
$\alpha$, stronger by tens of orders of magnitude than the thermal
photospheric emission. The other Lyman and Balmer line series also clearly
exceed the photospheric contribution. Notice also the clearly visible Lyman
and Balmer recombination continua. These are the smooth parts of the SED
between roughly 500 and 900 Å (with 921 Å the edge of our Lyman continuum at
$n>10$ instead of 912 Å for $n=\infty$) and 1000 and 4000 Å, respectively. The
continua for the other series are hidden by the photospheric contribution and
are thus negligible. In any case, as we discuss in Section 6.3, the strength
and shape of the continua are approximate—within a factor of a few, but also
possibly exceeding this—and are only meant to provide some guidance. Also,
note that the continua (hydrogen recombination and H-) are thought to mainly
come from the deeper and cooler region, the heated photosphere after (or
below) the end of our calculation.

###  3.2 Grid of models

We now present results from a large grid in accretion rate and mass,

$\displaystyle 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}\lesssim$
$\displaystyle~{}\dot{M}\lesssim 10^{-4}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$
(1) $\displaystyle 0.5~{}M_{\mathrm{J}}\lesssim$
$\displaystyle~{}M_{\mathrm{p}}\lesssim 20~{}M_{\mathrm{J}}.$ (2)

This wide parameter space is chosen to cover present and future observations,
which might reveal a population of closer-in accreting planets (Close, 2020).
The $R_{\mathrm{p}}(\dot{M},M_{\mathrm{p}})$ and
${T_{\textnormal{eff}}}(\dot{M},M_{\mathrm{p}})$ are given by the relations of
Section A.1, where ${T_{\textnormal{eff}}}$ is chosen so that the total flux
in the SED is equal to the sum of the internal and the incoming kinetic energy
flux (Equation (A5)). We take as the standard case the radius fit from the
“accretion hot-start” population and set $f_{\mathrm{fill}}=1$ for simplicity.

The upper accretion rate of $\dot{M}\approx
10^{-4}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ represents the highest values in
the runaway phase (Phase II) of classical in formation calculations, with a
dependence on the viscosity and the scale height of the PPD, as reviewed in
Helled et al. (2014). Included within this bound are thus the common maximum
values near $10^{-2}~{}M_{\mathrm{E}}\,\mathrm{yr}^{-1}\approx
10^{-4.5}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ (Bodenheimer et al., 2000;
Marley et al., 2007; Lissauer et al., 2009; Mordasini et al., 2012a; Tanigawa
& Tanaka, 2016; see also Schulik et al. 2020) so that our upper value
represents a conservative choice. It also equals the typical accretion rate
through the PPD $\dot{M}_{\mathrm{disk}}\sim
10^{-8}~{}{M_{\odot}}\,\mathrm{yr}^{-1}$ (Hartmann et al., 2016); the planet
would be intercepting the full typical $\dot{M}_{\mathrm{disk}}$ or a smaller
fraction of a higher global value.

The smallest $\dot{M}=10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ will turn
out to be roughly the lower limit needed to explain the PDS 70 b and c
observations of Haffert et al. (2019). Also, towards low $\dot{M}$ we expect
the photospheric noise to dominate over the accretion lines, in addition to
the total line luminosity becoming small, making it a less interesting case to
study.

The lower mass of $M_{\mathrm{p}}=0.5~{}M_{\mathrm{J}}$ corresponds to a free-
fall velocity near $30~{}\mathrm{km}\,\mathrm{s^{-1}}$ (see Figure 12), which
is the lower limit on $v_{0}$ for hydrogen line emission when the preshock gas
is in molecular form (Aoyama et al., 2018). As discussed in Section A.2.1, it
is not certain at these masses $\lesssim M_{\mathrm{J}}$ to what extent the
accreting gas is indeed in free-fall, and the preshock velocity could be
smaller. Thus this is an optimistic choice, especially since extinction by the
upper layers of the PPD could be important for small masses, which may not
carve out a deep gap.

Finally, we take as an upper limit $M_{\mathrm{p}}=20~{}M_{\mathrm{J}}$ for a
few reasons. One is to focus on objects that are predominantly formed as
planets, during the formation of which an accretion shock should occur,
whereas this is less clear for brown dwarfs (see discussion in Section 4.2 of
Baruteau et al. 2016). The planet and brown-dwarf mass functions
overlap222Core accretion can form objects up to several tens of
$M_{\mathrm{J}}$, with a low frequency (Mordasini et al., 2012a; Emsenhuber et
al., 2020b). and have a minimum near $M_{\mathrm{p}}=25~{}M_{\mathrm{J}}$
(Reggiani et al., 2016) before increasing towards low masses (Nielsen et al.,
2019). Thus most objects with $M_{\mathrm{p}}\lesssim 20~{}M_{\mathrm{J}}$ are
likely to have formed by core accretion (Schlaufman, 2018; Wagner et al.,
2019) and thus to have experienced an accretion shock, making them more
observationally relevant. Secondly, we will see that towards higher masses,
the line fluxes are relatively insensitive to the mass; thus stopping at 20 or
30 $M_{\mathrm{J}}$ makes little difference.

In any case, we emphasise that the range of $\dot{M}$ and $M_{\mathrm{p}}$ in
Equations (1) is not a prediction but rather meant as a conservative choice
for the input parameters, that is to say, a generous range of possibly
relevant values. We are not making statistical predictions for the H $\alpha$
luminosities as in Mordasini et al. (2017).

Figure 3:  SEDs of accreting planets for a grid of accretion rates and masses:
photospheric emission only (black lines) or with the shock-heated gas emission
(red). We took $f_{\mathrm{fill}}=1$, a distance of 150 pc, and the hot-start
fits (Section A.1.1) for the planetary radius $R_{\mathrm{p}}$ and effective
temperature ${T_{\textnormal{eff}}}$ (see subpanels), rounded to the nearest
100 K or up to ${T_{\textnormal{eff}}}=1200$ K. No contribution from a CPD
(Zhu, 2015) is shown. The ISM will absorb the Lyman lines and continuum but
locally they will affect the disc chemistry (Cleeves et al., 2015; Rab et al.,
2019).

Figure 3 shows a grid of SEDs for $M_{\mathrm{p}}/M_{\mathrm{J}}=1$, $3$,
$10$, $20$ and $\dot{M}/(M_{\mathrm{J}}\,\mathrm{yr}^{-1})=10^{-8}$,
$10^{-7}$, $10^{-6}$, $10^{-5}$. The peak intensities of the H $\alpha$ and
other hydrogen lines are significant relative to photospheric emission and
increase with $v_{0}$. In all panels except
$(\dot{M}=10^{-5}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1},M_{\mathrm{p}}=1~{}M_{\mathrm{J}})$,
the Ly $\alpha$ line at $\lambda=1215$ Å has the highest peak value. The
effective temperature ${T_{\textnormal{eff}}}$ is almost a monotonic function
of both $\dot{M}$ and $M_{\mathrm{p}}$ (cf. Section A.1.2). With increasing
planet mass and decreasing accretion rate, the Lyman continuum blueward of
$\lambda\approx 912$ Å becomes stronger relative to the hydrogen lines. As
discussed in Section 3.1.1, this means that a large fraction of the hydrogen
is ionized and that, in reality, a large fraction of the energy should be
converted into He and metal lines instead. However, again, this does not
affect the strength of the hydrogen lines.

In all cases shown, Lyman and Balmer lines have significant peaks above the
photospheric emission, because the peak of photospheric emission is at a
longer wavelength than these lines.

Figure 3 also shows that the ratio of the Ly $\alpha$ to the H $\alpha$ peaks
increases with planet mass, and that for $M_{\mathrm{p}}\approx
1~{}M_{\mathrm{J}}$, the ratio also increases with decreasing $\dot{M}$. This
is because at high postshock temperatures $T_{1}$ (high $M_{\mathrm{p}}$; see
Equation (A17)), hydrogen excitation occurs. This increases the abundance of
the absorber of H $\alpha$, namely electrons in the $n=2$ state, while
depopulating the absorber of Ly $\alpha$, electrons in $n=1$. This leads to a
lower H $\alpha$/Ly $\alpha$ ratio. Towards high postshock densities (high
$\dot{M}$; see Equation (A9)), both H $\alpha$ and Ly $\alpha$ are more
strongly absorbed. However, this absorption occurs in the upper regions (small
$\Delta z$), where the temperature is high but the excitation degree is low.
Normally, hotter gas emits and cooler gas absorbs, but since the hot gas has a
low excitation degree, the hot gas can absorb. This is a non-equilibrium
(NLTE) effect not captured by a time-independent approach. Therefore, since
the lower levels of hydrogen are more populated, Ly $\alpha$ absorption is
stronger than H $\alpha$ absorption. This leads to the increase of H
$\alpha$/Ly $\alpha$ with $\dot{M}$.

In Figure 3, we also see that longer-wavelength series (e.g., Paschen or
Brackett) are embedded in the photospheric signal but tend to emerge towards
larger masses and accretion rates. This suggests that high-resolution
spectroscopy of strongly accreting or massive planets might be able to detect
lines from these other series (see also Sections 5.3 and 5.4), unless infra-
red excess from dust particles in CPD is significant enough.

Particularly with the hot-start radii, the Mach number is not monotonically
proportional to $\dot{M}$ because of the non-monotonic dependence of $\dot{M}$
on the radius. Although increasing the mass flux of accreting gas $\dot{M}$
increases the amount of shock-heated (and thus emitting) gas, it is associated
with a larger planet radius at the same time. In turn, the larger planet
radius leads to a slower free-fall velocity at the planet surface $v_{0}$.
Figure 3 reflects this, given that, as we verified, the H $\alpha$ continuum
is an at least roughly monotonic function of ${T_{\textnormal{eff}}}$ at fixed
planet mass.

As one of the most important results of this work, Figure 4 shows the H
$\alpha$ line luminosity as a function of $\dot{M}$ and $M_{\mathrm{p}}$. We
extrapolated the model results down to $\dot{M}\sim
10^{-9.5}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$. Note that Figure 3 of Aoyama &
Ikoma (2019) is similar but was made for a fixed $R_{\mathrm{p}}$ = 2
$R_{\mathrm{J}}$, independently of $\dot{M}$ and $M_{\mathrm{p}}$. The H
$\alpha$ luminosity $L_{\mathrm{H}\,\alpha}$ ranges from $\sim
10^{-9}~{}{L_{\odot}}$ to $\sim 10^{-3}~{}{L_{\odot}}$ over the grid, overall
increasing monotonically with both $\dot{M}$ and $M_{\mathrm{p}}$. The
contours show that for $M_{\mathrm{p}}\gtrsim 3$–$5~{}M_{\mathrm{J}}$, the H
$\alpha$ luminosity is independent of $M_{\mathrm{p}}$ and is roughly linearly
proportional to $\dot{M}$. The first part of the reason for this is that the H
$\alpha$ luminosity turns out to be roughly linearly proportional to the
incoming kinetic energy flux
$L_{\mathrm{acc}}=GM_{\mathrm{p}}\dot{M}/R_{\mathrm{p}}$ (Aoyama et al.,
2018), especially at a fixed mass (see Figure 10 in Section A.1.2). The second
part is a simple one: the mass coordinate is on a linear scale, with only a
limited range ($\approx 1.3$ dex) relevant to planetary detections, while the
accretion rate axis is logarithmic and chosen to cover five orders of
magnitude.

Our $L_{\mathrm{H}\,\alpha}(\dot{M},M_{\mathrm{p}})$ relation (Figure 4) is
robust to changes in the model choices. We compare in Appendix D the
luminosity as obtained with the hot- and the cold-start $R_{\mathrm{p}}$ and
${T_{\textnormal{eff}}}$ relationships and find very little difference.
Similarly, varying $f_{\mathrm{fill}}$ from $f_{\mathrm{fill}}=1$ to
$f_{\mathrm{fill}}=0.01$ (not shown) changes the H $\alpha$ fluxes by at most
a factor of two333 At extremely low $f_{\mathrm{fill}}\lesssim 10^{-4}$, self-
absorption becomes very important. For more moderate values
$f_{\mathrm{fill}}\gtrsim 10^{-3}$ as inferred for young accreting stars
(Ingleby et al., 2013), self-absorption is not a significant effect.. This is
because the incoming gas mass at the shock ($\dot{M}$) is independent of
$f_{\mathrm{fill}}$, and H $\alpha$ emission is almost proportional to the
mechanical energy of incoming gas (Aoyama et al., 2018), as mentioned above.
See also the discussion in Section 6.2.

Figure 4 is meant as a tool for interpreting H $\alpha$ detections in terms of
fundamental planet parameters. Therefore, we also compare the luminosities
with those of a few low-mass objects (labeled contours), which we discuss in
the next section. Section 5.3 presents the Br $\gamma$, Pa $\alpha$, Pa
$\beta$, and H $\beta$ luminosities in a similar fashion.

Figure 4: Non-extincted H $\alpha$ luminosity from the planet-surface shock as
a function of accretion rate and planet mass (colorscale). Thin gray lines
highlight $\log L_{\mathrm{H}\,\alpha}/{L_{\odot}}=-9$ to $-4$ in steps of 1
dex. Shaded bands show non-dereddened 1-$\sigma$ contour regions for PDS 70 b
(dashed blue line: $10^{-5.9}~{}{L_{\odot}}$; Wagner et al., 2018; solid:
$10^{-6.8}~{}{L_{\odot}}$; Haffert et al., 2019) and PDS 70 c (green:
$10^{-7.1}~{}{L_{\odot}}$; Haffert et al., 2019). The contour for the
Hashimoto et al. (2020) value of
$L_{\mathrm{H}\,\alpha}=10^{-6.5}~{}{L_{\odot}}$ (not shown) would lie between
the two blue bands. The gray region is for the less secure protoplanet
candidate LkCa 15 b (black: $10^{-4.1}~{}{L_{\odot}}$; Sallum et al., 2015,
but see text).

##  4 Application to observations

In this section, we apply our results to a few (candidate) protoplanets444A
recent addition to this list may be 2M0249 c (catalog ) (Chinchilla et al.,
2021), with the caveat that chromospheric activity could be contributing to
the H $\alpha$.. Implications of the non-detections in dedicated surveys
(Cugno et al., 2019; Zurlo et al., 2020; Xie et al., 2020) are discussed in
Sections 6.5.

###  4.1 PDS 70 b and PDS 70 c

####  4.1.1 Brief partial review of observations

Wagner et al. (2018) reported the detection using MagAO of an H $\alpha$
signal from PDS 70 b, a companion in the gap in the transitional disc around a
young ($5.4\pm 1.0$ Myr; Müller et al. 2018) pre-main sequence
0.9–1.0-${M_{\odot}}$ star (Keppler et al., 2019; Wang et al., 2021)
discovered by Keppler et al. (2018). Then, the H $\alpha$ detection at PDS 70
b was confirmed by Haffert et al. (2019) using VLT/MUSE (Bacon et al., 2010).
They also reported the discovery in H $\alpha$ of PDS 70 c (catalog ), a
companion at the edge of the gap. From new VLT/SINFONI K-band data from
Christiaens et al. (2019b), Christiaens et al. (2019a) inferred the presence
of a circumplanetary disc around PDS 70 b, the first observational evidence
for a disc around a planet in a circumstellar disc. Using the near-infrared
(NIR) SED and the models of Eisner (2015), they derived an accretion rate
$\dot{M}\sim 10^{-7.5}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$. Also, Wang et al.
(2020) observed this system with Keck/NIRC2 and estimated the mass accretion
rate to be $\dot{M}=(3$–$8)\times 10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$
by comparing to the luminosity-evolution model of Ginzburg & Chiang (2019).
More recently, Stolker et al. (2020a) added the first detection of PDS 70 b at
4–5 $\upmu$m and re-analyzed the other data from 1 to 4 $\upmu$m, confirming
the finding by Wang et al. (2020) that a blackbody fits well the SED. Given
their modeling results, they concluded that PDS 70 b is likely surrounded by
some dusty material, which nevertheless lets (some) H $\alpha$ pass through.
Finally, thanks to $R\approx 500$ $K$-band spectra and astrometry of PDS 70 b
and c from VLTI/GRAVITY, Wang et al. (2021) found statistical support for a
more complex (non-blackbody) SED and for a small mass for PDS 70 b
($M_{\mathrm{p}}<10~{}M_{\mathrm{J}}$, but likely even lighter).

The H $\alpha$ signal of PDS 70 b has been detected with two different
instruments, with different luminosity determinations. From Wagner et al.
(2018), the luminosity can be derived as $L_{\mathrm{H}\,\alpha}=(1.4\pm
0.6)\times 10^{-6}~{}{L_{\odot}}$ following their data and description (see
Appendix E for details); this agrees with $L_{\mathrm{H}\,\alpha}=(1.3\pm
0.7)\times 10^{-6}~{}{L_{\odot}}$ derived by Thanathibodee et al. (2019). As
for Haffert et al. (2019), they obtain $L_{\mathrm{H}\,\alpha}=(1.6\pm
0.1)\times 10^{-7}~{}{L_{\odot}}$ for PDS 70 b and $(7.6\pm 1.3)\times
10^{-8}~{}{L_{\odot}}$ for PDS 70 c. Thus the luminosity for PDS 70 b derived
under the assumption of no extinction within the system, and ignoring the ISM
contribution of $A_{\mathrm{R}}\approx 0.02$–0.12 mag (see Appendix E), is
about lower by one order of magnitude than in Wagner et al. (2018).

Hashimoto et al. (2020) improved the data-correction method of the VLT/MUSE
data and estimated higher values $L_{\mathrm{H}\,\alpha}=(3.3\pm 0.1)\times
10^{-7}~{}{L_{\odot}}$ and $(1.3\pm 0.1)\times 10^{-7}~{}{L_{\odot}}$ for PDS
70 b and c, respectively. The value for PDS 70 b is still lower than in Wagner
et al. (2018) by a factor of four. This could be due to intrinsic variability
in the H $\alpha$ emission from PDS 70 b and/or from the known variability of
the star in the R band, combined with the way the contrast is measured.
However, Haffert et al. (in prep.) report that for a dozen MUSE measurements
over a period of three months, there is no variability in the H $\alpha$ flux
at the $\approx 30$% level. More recently, Zhou et al. (submitted) detected
PDS 70 b at H $\alpha$ with the Hubble Space Telescope (HST) and also found no
evidence for variability over an almost five-month baseline. Thus differences
in the data reduction seem to be a likely explanation for the differences.

####  4.1.2 No extinction at PDS 70 b and c?

While we assumed no extinction in this paper, Hashimoto et al. (2020) found
the H $\alpha$ from PDS 70 b and c to be likely strongly extincted
($A_{\mbox{H\,$\alpha$}}>2.0$ mag for PDS 70 b) based on the spectral width of
the H $\alpha$ line and their upper limit on the H $\beta$/H $\alpha$
fraction. Repeating their analysis with a more realistic opacity law, Marleau
et al. (subm.) derive even stronger constraints
($A_{\mbox{H\,$\alpha$}}\gtrsim 4$–8 mag for PDS 70 b).

However, if the observed spectral width is overestimated due to the
instrumental resolution (Thanathibodee et al., 2019), other solutions without
extinction are allowed (see Figure 3 in Hashimoto et al. 2020): towards lower
$n_{0}$ and $v_{0}$, both the flux ratio H $\beta$/H $\alpha$ and the line
widths are smaller, and there are matching combinations with
$A_{\mbox{H\,$\alpha$}}=0$ and smaller $f_{\mathrm{fill}}$. Thus, our
assumption can be consistent with the observational results. This solution
without (or with weak) extinction is preferred by the mass estimate of PDS 70
b and c (Bae et al., 2019; Stolker et al., 2020a; Wang et al., 2021). To
confirm whether the H $\alpha$ from PDS 70 b is strongly extincted, follow-up
observations with a higher spectral resolution are needed.

####  4.1.3 Derived accretion rate

Given these H $\alpha$ luminosities, our model yields
$M_{\mathrm{p}}$–$\dot{M}$ relations, shown in Figure 4 as line contours: blue
dashed (PDS 70 b; Wagner et al., 2018), blue solid (PDS 70 b; Haffert et al.,
2019), and green (PDS 70 c; Haffert et al., 2019), respectively. If
$M_{\mathrm{p}}=5$–$9~{}M_{\mathrm{J}}$ for PDS 70 b as (Wagner et al., 2018)
estimated and $M_{\mathrm{p}}=4$–$12~{}M_{\mathrm{J}}$ for PDS 70 c (Haffert
et al., 2019), our model implies $\dot{M}=(8.0\pm 4.8)\times 10^{-8}$ for PDS
70 b’s $L_{\mathrm{H}\,\alpha}$ from Wagner et al. (2018), $\dot{M}=(1.1\pm
0.3)\times 10^{-8}$ for PDS 70 b’s $L_{\mathrm{H}\,\alpha}$ from Haffert et
al. (2019), and $\dot{M}=(6.3\pm 3.1)\times
10^{-9}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ for PDS 70 c’s
$L_{\mathrm{H}\,\alpha}$ from (Haffert et al., 2019), respectively. If instead
$M_{\mathrm{p}}\lesssim 5~{}M_{\mathrm{J}}$, which is preferred for PDS 70 b
(Wang et al., 2021) and c (Bae et al., 2019) the constraints on $\dot{M}$ and
$M_{\mathrm{p}}$ are more accurately given in a joint form:
$\dot{M}M_{\mathrm{p}}\approx 5\times
10^{-7}~{}M_{\mathrm{J}}^{2}\,\mathrm{yr}^{-1}$ for the
$L_{\mathrm{H}\,\alpha}$ from Wagner et al. (2018), so that $\dot{M}\approx
1.7\times 10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ for
$M_{\mathrm{p}}=3~{}M_{\mathrm{J}}$, and, for the $L_{\mathrm{H}\,\alpha}$
from Haffert et al. (2019), $\dot{M}M_{\mathrm{p}}\approx 6.3\times
10^{-8}~{}M_{\mathrm{J}}^{2}\,\mathrm{yr}^{-1}$, implying $\dot{M}\approx
2.1\times 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ at
$M_{\mathrm{p}}=3~{}M_{\mathrm{J}}$. Towards low masses, these results depend
somewhat more on the choice of the radius fit (see Figure 15), but the main
source of uncertainty is the one in the observed value of
$L_{\mathrm{H}\,\alpha}$.

By extrapolating the empirical $L_{\mathrm{acc}}$–$L_{\mathrm{H}\,\alpha}$
relationship for Young Stellar Objects (YSOs) from Rigliaco et al. (2012),
Wagner et al. (2018) estimated $\dot{M}\approx
10^{-9}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ for PDS 70 b. Also, with the
$\dot{M}$–linewidth relationship of Natta et al. (2004), Haffert et al. (2019)
derived $\dot{M}\approx 2\times 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$.
Thus, applying stellar accretion models to planetary-mass observations yields
a lower mass accretion rate than from our model by a few orders of magnitude.
To estimate mass accretion rates, we suggest that our model constructed for
planet accretion should be used rather than the extrapolation of empirical
relationship from pre-main-sequence star studies. This is discussed briefly in
Section 6.1 but in more details in a companion publication, Aoyama et al.
(submitted).

Also, Thanathibodee et al. (2019) constructed a model of H $\alpha$ emission
focusing on PDS 70 b. They modeled the accreting gas as the source of the H
$\alpha$ rather than the postshock region that is the subject of this paper.
The accretion rate they estimate, $\dot{M}\approx~{}10^{-8.0\pm
0.6}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$, is larger than the results of
empirical $L_{\mathrm{H}\,\alpha}$–$L_{\mathrm{acc}}$ relationships and in
agreement with our results within the margin of error. As discussed in Section
6.1, whether the H $\alpha$ emission originates from the accretion flow or the
postshock region depends on whether the accreting gas is hot enough to emit H
$\alpha$. In fact, for PDS 70 b a contribution from both cannot be excluded
(Aoyama et al., submitted).

Finally, the upper limits on the Br $\alpha$ (Stolker et al., 2020a) and Br
$\gamma$ (Christiaens et al., 2019b; Wang et al., 2021) emission are discussed
in Section 5.4.

In summary, given an $L_{\mathrm{H}\,\alpha}$ measurement, our model yields
joint constraints on $\dot{M}$ and $M_{\mathrm{p}}$ at low masses, which seem
more likely for PDS 70 b and for PDS 70 c. (For higher masses,
$L_{\mathrm{H}\,\alpha}$ becomes relatively independent of $M_{\mathrm{p}}$.)
The uncertainty in $\dot{M}$ is $\approx 1$ dex. Our estimated $\dot{M}$
values for PDS 70 b and c are smaller than previously determined in the
stellar literature and similar to the results of Thanathibodee et al. (2019).
This is however a coincidence, since the two models have a very different
physical basis and predict in general distinct
$\dot{M}$–$L_{\mathrm{H}\,\alpha}$ relationships (Aoyama et al., submitted).
The main limitations on determining $\dot{M}$ are the uncertainties on the
line-integrated fluxes and the line widths, as well as the uncertainties about
the true masses.

###  4.2 LkCa 15 b

Following the discovery of a companion to LkCa 15 A (catalog ) by Kraus &
Ireland (2012), Sallum et al. (2015) reported the infrared detection of
further sources in the system using sparse-aperture masking (SAM).
Intriguingly, they also measured an H $\alpha$ signal which seemed to
originate at the position of LkCa 15 b. On the other hand, Thalmann et al.
(2016) analyzed scattered light from the disk and showed that the infrared
detections of the planetary candidates around LkCa 15 A could be false
positives related to features of the disc in scattered light. In addition,
observations by Mendigutía et al. (2018) using spectro-astrometry suggest that
the H $\alpha$ emission may not be coming from a point source but rather from
an extended region similar in size to the orbit of the claimed planet LkCa 15
b. Recently, Currie et al. (2019) conducted the first direct-imaging
observations of the LkCa 15 system. They provided evidence that there is no
point source at the location of the claimed planet (nor of the possible
further companions) but that in fact the SAM signal originates from disc
emission.

Despite the debate as to its origin, we will briefly analyze the H $\alpha$
signal at the position of a putative companion to LkCa 15 A as originating
from an accretion shock on the planet surface. The reported H $\alpha$
luminosity is $L_{\mathrm{H}\,\alpha}=10^{-4.1\pm 0.1}~{}{L_{\odot}}$ from
Sallum et al. (2015) but using the updated Gaia distance determination of
158.8 pc (Gaia Collaboration et al., 2018). From Figure 4 and assuming
$M_{\mathrm{p}}=10\,M_{\mathrm{J}}$, $\dot{M}=4.0^{+2.5}_{-0.1}\times
10^{-6}\,M_{\mathrm{J}}\,\mathrm{yr}^{-1}$. This accretion rate is not
implausible for a claimed forming gas giant, especially if it were undergoing
an accretion outburst.

Using instead the Rigliaco et al. (2012) approach as in Sallum et al. (2015)
and again with $M_{\mathrm{p}}=10~{}M_{\mathrm{J}}$ as an example, yields
$\dot{M}=3\times 10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ for
$R_{\mathrm{p}}=1.6~{}R_{\mathrm{J}}$ as Sallum et al. (2015) assumed. At this
$(\dot{M},M_{\mathrm{p}})$, our fits (Section A.1.1) yield
$R_{\mathrm{p}}=1.9~{}R_{\mathrm{J}}$ ($R_{\mathrm{p}}=1.4~{}R_{\mathrm{J}}$)
for the hot (cold) population, so that $1.6~{}R_{\mathrm{J}}$ is a reasonable
value, albeit perhaps on the small side. The upshot of the comparison is that
the $\dot{M}$ implied by the Rigliaco et al. (2012) relationship is one order
of magnitude smaller than derived with our approach; for PDS 70 b, it was a
few orders of magnitude. As discussed in Aoyama et al. (submitted), we suggest
that our models, which are tailored for the planetary case, should be used
instead of extrapolations from the stellar regime.

##  5 Further observational aspects

We now discuss to what extent high-mass and high-$\dot{M}$ planets can be
distinguished (Section 5.1) and the planet surface shock from the CPD shock
(Section 5.2), before presenting the line strengths and line ratios for
accretion-generated hydrogen lines other than H $\alpha$ (Section 5.3).
Finally, we discuss what information may be obtained from combining
observations of several lines for the same object (Section 5.4).

###  5.1 Distinguishing massive planets and strongly accreting planets from
the line profile?

Figure 5:  Spectral profile of H $\alpha$ for $v_{0}=100$, $150$,
$200~{}\mathrm{km}\,\mathrm{s^{-1}}$ at $n_{0}=10^{11}$ cm-3 (left panel) and
for $n_{0}=10^{9}$, $10^{11}$, $10^{13}~{}\mathrm{cm}^{-3}$ at
$v_{0}=150~{}\mathrm{km}\,\mathrm{s^{-1}}$ (right panel), typical for
planetary masses (see Figure 12). Each line is normalized by its peak value.
The top axis shows the Doppler shift velocity, even though some features come
from natural and not Doppler broadening. The spectral resolution of MUSE at H
$\alpha$ is indicated ($R=2516$; red bar). The H $\alpha$ line width depends
on $n_{0}$ more than on $v_{0}$ because high $n_{0}$ leads to self-absorption
near the peak.

When characterizing gas giants from their H $\alpha$ luminosity, their mass
and mass accretion rate are degenerate because the luminosity depends on their
product. However, this degeneracy can be lifted by spectroscopic observations
of H $\alpha$, which was demonstrated by Aoyama & Ikoma (2019) in the case of
PDS 70 b and c.

Recall that the preshock velocity $v_{0}$ mainly depends on $M_{\mathrm{p}}$,
while the number density $n_{0}$ is mainly set by $\dot{M}$ (see Figure 13).
Figure 5 shows H $\alpha$ line shapes for several values of $v_{0}$ and
$n_{0}$. The preshock velocity $v_{0}$ sets the shock strength—and thus the
temperature just after the shock (Equation (A17))—but barely the line width,
which is mainly set by Doppler broadening. This is because the gas that
becomes ionized and then recombines is the main source of H $\alpha$, and not
the gas immediately after shock, even if there is some amount of excited
hydrogen there (see down to $\Delta t\approx 10^{-3}$ s in Figure 1). Note
that Figure 5a a small redshift (of several $\mathrm{km}\,\mathrm{s^{-1}}$)
can be seen in the line profiles. This is due to the non-zero settling speed
of the emitting gas.

The line profile can be divided into three parts: a narrower Gaussian, a
broader Gaussian, and a Lorentzian profiles, with the latter visible further
from the line center. The layers that emit the two Doppler profiles are
separated by the highly ionized region (at $\Delta t\approx 10^{-3}$–10 s in
Figure 1). The H $\alpha$ intensity coming from the deeper layers is much
larger. Consequently, even the width of the line where the energy density is
10 % of the maximum ($W_{10}$ in e.g. Thanathibodee et al. 2019) reflects only
the narrower Doppler component, coming from the hydrogen-ion recombination
region at low temperatures, in the deep layers. The gas temperature at which
hydrogen recombination begins barely depends on $v_{0}$, because it
corresponds to the hydrogen ionization energy of $13.6$ eV. We can see the
thermal broadening of the optically thinner gas just after the shock only far
from the line center, at $\lambda\gtrsim 6564$ Å and $\lesssim 6562$ Å (i.e.,
$|\Delta v|\gtrsim 50~{}\mathrm{km}\,\mathrm{s^{-1}}$ away from the shock).
Since the hot gas immediately below the shock has a high velocity and is
travelling away from the observer, the red half of the line is more broadened
by this mechanism than the blue half. However, as seen in the right panel, the
resolution of MUSE ($R=2516$ at H $\alpha$; Eriksson et al., 2020) is not
sufficient to distinguish this.

As shown in the right panel of Figure 5, $n_{0}$ changes the width of the
normalized line dramatically. However, increasing $n_{0}$ hardly broadens the
H $\alpha$ line because the pressure broadening is negligible relative to the
natural broadening. A high $n_{0}$ leads to H $\alpha$ self-absorption in the
postshock gas (in the top part of the flow), which flattens the line peak.
Since we normalized the line flux at the peak, the self-absorbed line looks
broader (see Aoyama et al., 2018 for the non-normalized profile). However, the
lines for higher $n_{0}$ are brighter than lower ones despite the absorption.
This effect becomes significant for $n_{0}>10^{11}$ cm-3 in the right panel.
Also, a lower $n_{0}$ ($10^{9}$ cm-3) can lead to slightly broader profile.
While the wider component that comes from the shallower region is independent
of the density, the narrower component that comes from the deeper region gets
weaker with decreasing $n_{0}$ due to a lower excitation degree. Thus, the
wider component gets stronger relative to the narrower component.

As shown in Figure 5, the current spectral resolution of MUSE is not enough to
distinguish the profiles clearly, while it barely resolves the spectral
profile for higher density (Eriksson et al., 2020). However, it is not
possible to determine in general what minimum spectral resolution is required
for distinguishing high accretion rates from high masses because it depends on
the relative uncertainty in the flux as well as on the planet properties
through the dependence of the line profile on $(n_{0},v_{0})$.

###  5.2 Distinguishing planetary-surface and CPD-surface shocks?

Figure 6: SEDs for two different models of shock-heated gas emission:
planetary photospheric emission plus planetary-surface shock (red) and CPD
surface-shock (blue). The photosphere SED in the surface-shock case is also
shown for reference (black). The $\dot{M}$, $M_{\mathrm{p}}$, $R_{\mathrm{p}}$
values (see figure and text) are chosen such that
$L_{\mathrm{H}\,\alpha}\approx 1.5\times 10^{-6}~{}{L_{\odot}}$ in both cases.
The energy flux is calculated for the distance of 150 pc. Left panel: Global
SED. Right panel: H $\alpha$ line profile against Doppler shift velocity from
the line center. The photospheric emission is not visible.

Hydrodynamic simulations report that gas accreting toward proto-gas-giants
goes through multiple shocks (e.g. Kley, 1999; Tanigawa et al., 2012). The gas
that falls onto the CPD yields a shock, which can emit H $\alpha$ near the
planet. However, far from the planet, the shock is not strong enough to emit H
$\alpha$, and only the part of the shock close to the planet can contribute to
the H $\alpha$ emission (Aoyama et al., 2018). If the gas ultimately joining
the planet passes firstly through a shock at the surface of the CPD and
secondly through the planetary surface shock, the former is negligible for H
$\alpha$ emission, because most of the gas hits the CPD at the far region.
There, the free-fall velocity is too small for significant H $\alpha$
emission. For example, when the CPD is truncated near the planet (and the gas
is accreting by magnetospheric accretion or it is falling directly onto the
planet from the PPD), the planetary surface shock dominates the emission.
However, when most gas passes through boundary-layer accretion (see, e.g.,
Dong et al., 2020 rather than a planetary surface shock, the CPD surface shock
becomes significant. Thus it would be desirable to distinguish the source of
shock excess.

In Figure 6, we compare the SEDs from the planetary surface shock and the CPD
surface shock. The shock excess from the latter highly depends on the gas
accretion model. Here, guided by the results of an isothermal 3D hydrodynamic
simulation (Tanigawa et al., 2012), we assume the following:

  1. 1.

All gas accretes vertically from the protoplanetary disk onto the CPD with a
free-fall velocity starting from infinity set by the protoplanet’s gravity,
${v_{\textnormal{ff}}}(r)=\sqrt{2GM_{\mathrm{p}}/r}$, where $r$ is the radial
distance from the protoplanet’s center.

  2. 2.

The mass accretion flux onto the CPD is spatially constant within 0.1 Hill
radius, and zero outside of this. The inner disc radius does not matter
because the contribution of the outer region to the intensity dominates over
that of the inner region (Aoyama et al., 2018).

Non-isothermal simulations might find a different flow pattern (especially
concerning assumption 2.) but this should capture qualitatively the main
differences between the planet-surface and CPD-surface shock cases.

As an example, we consider system parameters that could be appropriate for PDS
70 b. The semi-major axis is $22.6$ au and the central-star mass is
$0.82~{}{M_{\odot}}$ (Riaud et al., 2006). The mass of PDS 70 b is very
uncertain, with some indications for a low mass of a few $M_{\mathrm{J}}$
(Mesa et al., 2019; Stolker et al., 2020a). Here, we choose a somewhat high
value within the range considered in the literature, namely
$M_{\mathrm{p}}=8~{}M_{\mathrm{J}}$ (Wagner et al., 2018). This places us in
the flat part of the $L_{\mathrm{H}\,\alpha}(\dot{M},M_{\mathrm{p}})$
contours, making the choice of $\dot{M}$ simple for the surface-shock case,
with $\dot{M}=10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$. To have the same H
$\alpha$ luminosity $L_{\mathrm{H}\,\alpha}\approx 1.5\times
10^{-6}~{}{L_{\odot}}$ in the CPD-shock case (Wagner et al., 2018), we take
for the CPD-shock case $\dot{M}=10^{-3}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$,
with the same mass. In both cases, we set
$R_{\mathrm{p}}=1.7~{}R_{\mathrm{J}}$ and ${T_{\textnormal{eff}}}=1600$ K
(Section A.1). For the surface-shock case, $f_{\mathrm{fill}}=1$ is also
assumed.

What fraction of the accreting gas in the CPD-shock case can produce hydrogen
lines? The chosen parameters yield ${R_{\textnormal{Hill}}}\approx
7000~{}R_{\mathrm{J}}=3.3$ au. From Figure 12, at this mass only the gas
within $r\approx 30~{}R_{\mathrm{J}}$ has ${v_{\textnormal{ff}}}\gtrsim
30~{}\mathrm{km}\,\mathrm{s^{-1}}$ and thus contributes to H $\alpha$. This
region corresponds to $0.005{R_{\textnormal{Hill}}}$, and thus 0.2 % of the
area over which the gas is assumed to accrete. Thus only this small fraction
of the accretion rate is available for producing H $\alpha$. This partly
explains the need for a high total $\dot{M}$ compared to the surface-shock
case to have the same $L_{\mathrm{H}\,\alpha}$. The other reason is that
$v_{0}$ is smaller everywhere on the CPD than on the planet surface (see
Figure 12), so that $\dot{M}$ must be higher to compensate since
$F_{\mathrm{H}\,\alpha}$, the H $\alpha$ flux at the object’s surface, roughly
scales with the incoming kinetic energy flux $F_{\mathrm{acc}}\propto
n_{0}v_{0}^{3}$.

The left panel of Figure 6 displays the global SEDs. The red and blue lines
correspond to the SED assuming planetary surface shock and CPD surface shock,
respectively. The black line corresponds to the SED without a shock excess for
reference. In the Lyman and Balmer continua, namely for $\lambda\lesssim 3000$
Å, the two SEDs differ by more than a factor of ten. This comes from the
difference in gas temperature after the two shocks. For the CPD surface shock,
the regions far from the planet dominate the shock excess emission because of
their large emitting area. This is a relatively weak shock, associated with a
low temperature.

This temperature difference can be also seen in the H $\alpha$ profile in
Figure 6b: the profile is narrower in the CPD surface-shock case than in the
planet surface-shock case when taking the different heights into account
(i.e., looking at the full width at half-maximum (FWHM)). The difference is
small but might be detectable in the future with high-resolution observations.
For example, for the case of Figure 6, the wavelength difference at the 10% of
the peak corresponds to $17~{}\mathrm{km}\,\mathrm{s^{-1}}$, which is
distinguishable with the spectral resolution of $R\approx 18,000$ (which is
much higher than the resolution of MUSE with $R\approx 2500$). At lower flux
levels relative to the peak, the difference is greater (e.g., $\approx
25~{}\mathrm{km}\,\mathrm{s^{-1}}$ at 5%) but it is not clear whether such
levels can be robustly extracted because of the high contrast relative to the
peak. The continuum emission from the CPD (Zhu, 2015) or from the planet’s
photosphere are likely not important, but the star’s chromosphere could
contribute (Manara et al., 2013, 2017; Venuti et al., 2019). The stellar
accretion-induced H $\alpha$ is possibly Doppler-shifted away from the
planet’s signal as in Haffert et al. (2019). However, the most important
limitation is probably the maximal contrast allowed by the instrument.

The narrower H $\alpha$ and the weaker recombination continua in the CPD case
means a weaker shock, which can also occur at the surface of a less massive
planet. However, in such a case, the density should higher than the CPD case,
and one can distinguish these two cases.

Our model does not include some continuum sources such as a heated photosphere
(e.g. Königl, 1991; Calvet & Gullbring, 1998) or, if present, the boundary
layer (e.g. Kenyon & Hartmann, 1987; Dong et al., 2020), which are well
modeled in the stellar accretion context. Such continuum sources can change
the spectral appearance, but it is unfortunately difficult to say how
important this would be.

In summary, for a given H $\alpha$ luminosity, the resolution of MUSE is not
sufficient to distinguish an accretion shock on the planetary surface from the
one on the CPD, but high-resolution spectroscopy might be able to do so.

###  5.3 Predictions for hydrogen lines other than H $\alpha$

The H $\alpha$ line on which we have focused so far is only one of the 55
hydrogen lines we model. Recently, Eriksson et al. (2020) reported an H
$\beta$ flux for the $\approx 10$-$M_{\mathrm{J}}$ companion Delorme 1 (AB)b
(catalog ) with the MUSE instrument on the VLT. Also, the upcoming, first-
light HARMONI555See https://harmoni-elt.physics.ox.ac.uk. integral field unit
(IFU; Thatte et al., 2016; Rodrigues et al., 2018) on the ELT is expected to
provide $R\approx 17,000$ spectroscopy between 0.8 and 2 $\mu$m (thus
including for instance Pa $\beta$ and Br $\gamma$); the second-generation
instrument HIRES for the ELT (Marconi et al., 2016, 2018; Tozzi et al., 2018)
will cover 1–1.8 $\upmu$m and thus should observe Paschen lines with the high
spectral resolution of $R\approx 100,000$; and the University of Tokyo Atacama
Observatory (TAO) should be able to detect Pa $\alpha$ thanks to its location
at 5,640 m (Yoshii et al., 2010). Finally, the Keck Planet Imager and
Characterizer (KPIC) (Jovanovic et al., 2019; Morris et al., 2020) aims at
obtaining $R=35,000$ spectroscopy in the $K$, $L$, and $M$ bands ($\approx
2$–5 $\upmu$m). Clearly, it is timely to extend the luminosity predictions to
lines other than H $\alpha$.

Figure 7: As in Figure 4 but for H $\beta$, Pa $\alpha$, Pa $\beta$, Pa
$\gamma$, Br $\alpha$, and Br $\gamma$. The color map and the solid contour
lines show the results from the radius fit to the warm population, and the
dashed contour lines use the cold population. Black lines show 5-$\sigma$
upper limits of $\log({L_{\textnormal{H}\,\beta}}/{L_{\odot}})<-7.03$
(Hashimoto et al., 2020) and
$\log({L_{\textnormal{Br}\,\gamma}}/{L_{\odot}})<-7.48$ (Christiaens et al.,
2019b) for PDS 70 b, and of
$\log({L_{\textnormal{Pa}\,\beta}}/{L_{\odot}})<-9.20$ for a putative planet
at 25 au in TW Hya (Uyama et al., 2017). Wang et al. (2021) set Br $\gamma$
upper limits for PDS 70 b and c similar to the line of Christiaens et al.
(2019b). In the Br $\alpha$ panel, black contours are for the detected flux at
NACO/NB4.05, $\log({L_{\textnormal{Br}\,\alpha}}/{L_{\odot}})=-5.88$ (Stolker
et al., 2020a), which is an upper limit on the line flux if the continuum is
much stronger, as is expected.

In Figure 7, we show the line luminosity of H $\beta$, Pa $\alpha$, Pa
$\beta$, and Br $\gamma$ for the same grid of $(\dot{M},M_{\mathrm{p}})$ as in
Figure 4. The luminosities range from $L_{\mathrm{line}}\sim 10^{-12}$ to
$10^{-4}~{}{L_{\odot}}$, increasing with $\dot{M}$ and $M_{\mathrm{p}}$ as for
H $\alpha$, and with the same qualitative shape of a very weak mass dependence
for $M_{\mathrm{p}}\gtrsim 3$–$5~{}M_{\mathrm{J}}$. The contours are very
similar using the fit to the hot- or cold-start populations. We show upper
limits for a few objects but discuss them below in Section 5.4.

Figure 8: Flux of H $\beta$, Pa $\alpha$, Pa $\beta$, Pa $\gamma$, Br $\alpha$
and Br $\gamma$ relative to H $\alpha$ as a function of the H $\alpha$ flux at
the surface of the object $F_{\mathrm{H}\,\alpha}$ and the preshock number
density $n_{0}=10^{9}$–$10^{14}$ cm-3 (colour). For each $n_{0}$, the preshock
velocity $v_{0}$ is varied from 20 to 200 $\mathrm{km}\,\mathrm{s^{-1}}$. The
H $\alpha$ luminosity is related to the $x$-axis by
$L_{\mathrm{H}\,\alpha}=4\times
10^{-8}F_{\mathrm{H}\,\alpha}(R_{\mathrm{p}}/1.5~{}R_{\mathrm{J}})^{2}(f_{\mathrm{fill}}/1)~{}{L_{\odot}}$.
The apparent flux $\mathfrak{F}_{\mathrm{H}\,\alpha}$ of an
1.5-$R_{\mathrm{J}}$ object at 100 pc is given by
$\log\mathfrak{F}_{\mathrm{H}\,\alpha}=\log F_{\mathrm{H}\,\alpha}-20.9$ dex.

Next, Figure 8 shows the intensity ratio of H $\beta$, Pa $\alpha$, Pa
$\beta$, and Br $\gamma$ relative to H $\alpha$ as a function of the H
$\alpha$ flux at the surface of the object. (Absolute fluxes for H $\beta$, Pa
$\alpha$, and Pa $\beta$ as a function of $(n_{0},v_{0})$ can be found in
Aoyama et al. (2018).) The ratios span H $\beta$/H $\alpha$ $\approx 0.1$–2 to
Br $\gamma$/H $\alpha$ $\approx 0.001$–0.05, with typically $X/$H $\alpha$
$\approx$0.03–0.3 for $X=$ H $\beta$, Pa $\alpha$, Pa $\beta$.

All four ratios are more or less flat within 0.5 dex at low H $\alpha$ flux
(within 0.7 dex for Br $\gamma$/H $\alpha$) but start to increase around an H
$\alpha$ flux $F_{\mathrm{H}\,\alpha}\sim
10^{7}\,\mathrm{erg\,s^{-1}\,cm^{-2}}$. For Pa $\alpha$ the rise is only
moderate. This change of slope occurs because the H $\alpha$ saturates due to
self-absorption in the shock-heated gas as $v_{0}$ increases, while the other
lines do not saturate. Thus towards high $F_{\mathrm{H}\,\alpha}$, $v_{0}$
must increase faster along the $x$ axis, which leads to stronger other lines
(and thus ratios) since they hardly display self-absorption. Also, when the
starting level of the transition is high (i.e., Pa $\beta$ or Br $\gamma$:
$n=5\rightarrow 3$ or $n=7\rightarrow 4$, respectively), some features at
weaker H $\alpha$ ($F_{\mathrm{H}\,\alpha}\lesssim
10^{6}~{}\mathrm{erg\,s^{-1}\,cm^{-2}}$) become more clearly visible. To
understand this, one should recall that for a given $n_{0}$, a low H $\alpha$
flux means a low gas temperature, and that for low temperatures,
$F_{\mathrm{H}\,\alpha}$ increases faster than the other lines due to its
lower excitation energy, and conversely at higher temperatures. Therefore, the
ratio first decreases and then increases towards high temperature
($F_{\mathrm{H}\,\alpha}$).

Also, our model predicts that the Ly $\alpha$ is much stronger than all lines
and carries most of the incoming shock energy, which has been converted to
radiation (see Figures 1 and 3). This is an important input for models of CPD
chemistry (Cleeves et al., 2015; Rab et al., 2019). However, interstellar
extinction is too strong for Ly $\alpha$ to be detected, excepting a few young
stellar systems such as TW Hydra. However, other lines induced by planetary Ly
$\alpha$ might be detected as for some accreting stars, including NIR
fluorescent molecular hydrogen lines (e.g., Herczeg et al., 2004).

Finally, Zhou et al. (submitted) detected for the first time PDS 70 b in the
Balmer continuum with the F336W ($U$-band) filter of the Wide-Field Camera 3
(WFC3) onboard HST. While our model does predict hydrogen recombination
continua, it is constructed with a focus on the hydrogen lines. The predicted
continua are less reliable and likely overestimated especially for massive
planets. This is because where the gas is almost ionized, the only effective
coolant in our model is the recombination continua, while in reality metal
lines should play an important role. Including them would decrease the fluxes
currently predicted in the continua. To treat them more accurately, a model
update is planned in the future.

###  5.4 Combining detections of multiple accretion lines

We briefly discuss the few current application of these predictions in the
planetary-mass regime. Firstly, the 5-$\sigma$ Br $\gamma$ upper limit for PDS
70 b from Christiaens et al. (2019b),
$\log{L_{\textnormal{Br}\,\gamma}}/{L_{\odot}}<-7.48$, is shown in Figure 7.
From $R\approx 500$ VLTI observations, Wang et al. (2021) confirmed this upper
limit to within 0.01 dex, and also placed a similar 5-$\sigma$ limit of
$\log{L_{\textnormal{Br}\,\gamma}}/{L_{\odot}}<-7.58$ for PDS 70 c. Taken at
face value, this implies for PDS 70 b that $\dot{M}\leqslant 2\times
10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ (and somewhat less for PDS 70 c)
if $M_{\mathrm{p}}\gtrsim 8~{}M_{\mathrm{J}}$, but is much less constraining
at lower $M_{\mathrm{p}}$.

Also, Stolker et al. (2020a) detected PDS 70 b for the first time in the Br
$\alpha$ filter of NACO/VLT (NB4.05; the effective width is 0.0616 $\upmu$m
and not 0.02 $\upmu$m666The too-narrow width came from
http://www.eso.org/sci/facilities/paranal/decommissioned/naco/inst/filters.html
(corrected in January 2021), which however provides the correct filter curve,
and is often quoted in the literature (e.g. Janson et al., 2008; Quanz et al.,
2010; Meshkat et al., 2014; Kervella et al., 2014; Stolker et al., 2020b, a).
It however does not change the results in those studies. The correct value is
from the SVO at http://svo2.cab.inta-
csic.es/theory/fps3/index.php?id=Paranal/NACO.NB405.). The Br $\alpha$ line
can be used as an accretion tracer (Komarova & Fischer, 2020), but in this
case the flux is consistent with the blackbody emission matching the global
SED. Namely, given the observed H $\alpha$ flux, the Br $\alpha$ is expected
to be embedded in the continuum (see Figure 8), and the upper limit of the Br
$\alpha$ emission from the shock itself is the observed line-integrated
luminosity of $L=1.3\times 10^{-6}~{}{L_{\odot}}$. This implies that
$\dot{M}\leqslant 10^{-5}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ if
$M_{\mathrm{p}}\gtrsim 5~{}M_{\mathrm{J}}$, which is a mass range not favoured
by the analysis of Wang et al. (2021) but not completely unlikely. For
$M_{\mathrm{p}}$ of at most a few $M_{\mathrm{J}}$, which seems more plausible
(Stolker et al., 2020a; Wang et al., 2021), the
${L_{\textnormal{Br}\,\alpha}}$ upper limit is effectively not constraining,
with $\dot{M}\leqslant 10^{-4}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$. Such high
rates are barely expected from the theoretical side. Even the somewhat lower
upper limits at high masses are consistent with the results from Figure 4,
$\dot{M}\approx 8\times 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ from the
Wagner et al. (2018) measurement or $\dot{M}\approx 1\times
10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ from Haffert et al. (2019).

For comparison, Hashimoto et al. (2020) derived from their re-analysis of
archival MUSE data an H $\beta$ flux upper limit of $2.3\times 10^{-16}$ erg
s-1 cm-2 for PDS 70 b, corresponding to ${L_{\textnormal{H}\,\beta}}<9.2\times
10^{-8}~{}{L_{\odot}}$ (3$\sigma$). This H $\beta$ is inconsistent with the H
$\alpha$ result (see Figure 4). The same is true for PDS 70 c. From this,
Hashimoto et al. (2020) concluded that there must be differential extinction.
As justified in Section 2.4, we emphasise that our results are without taking
extinction into account, which we do separately in Marleau et al. (subm.).

In the same class as CT Cha B/b (catalog ) or DH Tau B/b (catalog ), some of
the only low-mass putative accretors for which other accretion lines have been
detected (Schmidt et al., 2008; Bonnefoy et al., 2014; Zhou et al., 2014;
Bowler et al., 2014; Wu et al., 2015), is Delorme 1 (AB)b, for which Eriksson
et al. (2020) measured H $\beta$ (and also He i lines and upper limits on the
infrared Ca ii triplet). They infer a mass $M_{\mathrm{p}}\approx
12~{}M_{\mathrm{J}}$ and radius $R_{\mathrm{p}}=1.6~{}R_{\mathrm{J}}$ from
(hot-start) evolutionary models by combining their results to the photometry
of Delorme et al. (2013). They report
$\log(L_{\mathrm{H}\,\alpha}/{L_{\odot}})=-7.05\pm 0.06$, from which they
infer $\dot{M}\approx(0.8$–$3.0)\times
10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ by combining their
$L_{\mathrm{H}\,\alpha}$ measurement with different models. Combined with ours
(Aoyama et al., 2018; Aoyama & Ikoma, 2019; this work), they derive
$\dot{M}=1\times 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ and
$M_{\mathrm{p}}=11~{}M_{\mathrm{J}}$, using the radius
$R_{\mathrm{p}}=1.6~{}R_{\mathrm{J}}$ suggested by the photometry. From Figure
7, the predicted H $\beta$ luminosity is ${L_{\textnormal{H}\,\beta}}\approx
4\times 10^{-8}~{}{L_{\odot}}$. However, Eriksson et al. (2020) measured
${L_{\textnormal{H}\,\beta}}=10^{-8}~{}{L_{\odot}}$, which is a factor four
lower than the model prediction. This suggests that as for PDS 70 b,
extinction is affecting the measurement of accretion tracers. The interesting
difference is that there is up to now no evidence for an accretion disc around
Delorme 1 (AB)b, but the constraints on a disc mass or surface density are not
clear (Eriksson et al., 2020). Thus a detection of (or upper limit on) a disc
around Delorme 1 (AB)b, as well as more observational information on other
hydrogen lines would be useful. On the theoretical side, predictions for the
other currently available lines (He i and the infrared Ca ii triplet) would be
welcome.

Finally, in Figure 7, we show the 5-$\sigma$ upper limit on Pa $\beta$
emission for the TW Hya (catalog ) disc in its gap at 25 au,
$\log({L_{\textnormal{Pa}\,\beta}}/{L_{\odot}})<-9.20$, derived by Uyama et
al. (2017) using Keck. At 95 au, where there is an other gap, the upper limit
is $\log({L_{\textnormal{Pa}\,\beta}}/{L_{\odot}})<-9.79$ (not shown). This is
consistent with the mass constraints of $M_{\mathrm{p}}\lesssim
0.5~{}M_{\mathrm{J}}$ from van Boekel et al. (2017) for TW Hya. Looking at
Figure 7, the interesting implication is that Keck at Pa $\beta$ is sensitive
to planets with a relatively low mass or accretion rate, at least for the
nearest protoplanetary discs.

##  6 Discussion

We now take a critical look at different aspects of the model and results
presented in this work (Sections 6.1–6.4), and discuss the non-detections of
recent H $\alpha$ surveys (Section 6.5). Some caveats about the model were
already discussed in Section 4.2 of Aoyama et al. (2018) and we do not repeat
them here. We compare with other recent models of H $\alpha$ emission from
accreting planets (Thanathibodee et al., 2019; Szulágyi & Ercolano, 2020) and
discuss the validity of their approach in a different work (Aoyama et al.,
submitted). Appendix B already details why the physical assumptions behind
Storey & Hummer (1995) do not apply to the planetary surface shock.

###  6.1 Emission by the gas accreting onto the planet

Our model treats only the shock-heated gas. Several previous studies focusing
on stellar-mass objects considered the gas flowing onto the accretor as the
source of observed H $\alpha$ excess, assuming an unknown source of heating
there (e.g. Hartmann et al., 1994; Edwards et al., 2013), as opposed to
heating provided by the shock. Since the shock on accreting stars makes the
gas too hot ($>10^{6}$ K) to emit hydrogen lines, the shock-heated gas is
negligible for hydrogen-line emission. On the other hand, the planet-surface
shock emits significant hydrogen lines, while the accreting gas should be
cooler and emit weaker or no hydrogen lines compared to the stellar cases (see
also the discussion in Aoyama et al. submitted). This is the reason why we
neglect hydrogen line emission from the region.

Even if the gas in the accretion flow is too cool to emit H $\alpha$, the warm
gas plays a significant role in excess emission other than from hydrogen
lines. In the T Tauri-star context, Calvet & Gullbring (1998) found that a
Balmer recombination continuum is emitted by the gas accreting onto the object
rather than the postshock gas. The hydrogen recombination continua are well
modeled and compared with observational results in the stellar accretion
context. As seen in Figures 2 and 3, the planetary continua are predicted to
be several (up to tens of) orders of magnitude stronger than the contribution
from the planetary photosphere. Therefore, detecting in the planetary context
a hydrogen continuum much stronger than the photospheric emission would lend
support to our emission models.

###  6.2 Effect of the accretion geometry:
spherical versus magnetospheric

The predicted H $\alpha$ luminosity (see Figure 4) was derived explicitly in
the context of spherical accretion onto the protoplanet’s surface. However, it
also represents the signal expected for any accreting planet, regardless of
the accretion geometry, i.e., spherical or magnetospherical. In the case of
magnetospheric accretion with a filling factor $f_{\mathrm{fill}}<1$, relative
to spherical accretion, the kinetic-energy flux locally is higher by
$1/f_{\mathrm{fill}}$ but the accreting area is smaller by
$f_{\mathrm{fill}}$. These effects cancel each other out to a large extent
over most of parameter space. However, more precisely, spherical accretion
yields an upper limit to the H $\alpha$ intensity because of H $\alpha$ self-
absorption in the shock-heated gas. To emit more intense H $\alpha$ by
avoiding self-absorption, the gas needs to be less dense, which, at a given
$\dot{M}$, will be the case for higher $f_{\mathrm{fill}}$. Note that in a
realistic situation, the infalling gas and dust could absorb a part of the
flux emitted at the planet surface. This is explored systematically in Marleau
et al. (subm.).

Also, in the scenario of accretion onto a circumplanetary disk, the H $\alpha$
emission is at most roughly 1 % of the H $\alpha$ coming from the planet
surface for a similar planet mass, at least for the (simple) disc model and
scaling assumed in Aoyama et al. (2018). (In Section 5.2 we had compared the
line shapes at fixed total luminosity.) Therefore, when a strong shock occurs
on the planetary surface, regardless of the geometry, the CPD surface shock is
negligible.

###  6.3 Helium and metal lines

For accreting stars, helium and metal lines (He i, Ca ii, Na i, O i, etc.) are
also detected and used as indicators of stellar accretion as for hydrogen
lines (e.g. Kastner et al., 2002). At the upper edge of the planetary-mass
range, $M_{\mathrm{p}}\approx 20~{}M_{\mathrm{J}}$ leads to
${v_{\textnormal{ff}}}\approx 200$ km/s for $R_{\mathrm{p}}=2R_{\mathrm{J}}$
(see Figure 13), so that the postshock gas temperature can exceed
$T_{1}\approx 10^{6}$ K (Equation (A17)). In that case, metal lines instead of
hydrogen lines are responsible for the dominant emission processes at
ultraviolet (UV) wavelengths.

Our estimate of the hydrogen recombination continua would change when
including metal lines. When hydrogen ionization proceeds and neutral hydrogen
is minor (e.g. at $T\approx 10^{5}$ K in Figure 1), the gas should cool
through hydrogen recombination continua and/or metal lines. Presently, because
we do not include metal lines, almost all the thermal energy is converted into
the continua, so that they are overestimated. The hydrogen recombination
continua are mostly used as accretion indicator of protostars (e.g. Calvet &
Gullbring, 1998). For planetary accretion, they are currently expected to
highly exceed photospheric emission at UV wavelengths, as in the example in
Figure 2, but this should be re-assessed once more complete models are
available.

###  6.4 On the physical size of the cooling region

Towards high $v_{0}$, the line-forming region can be at a depth of order of
the planet size, which is clearly unrealistic. By contrast, in Figure 1 it is
at $\Delta z\approx 3\times 10^{7}$ cm, which represents $0.002R_{\mathrm{p}}$
and is much smaller than the planet size and thus reasonable. The reason for
the large extent of the cooling zone in some of the other cases is that there
is no cooling by helium nor atomic metals. In any case, we remind that we
always resolve the line-forming region because of adaptive time-stepping,
which keeps the relative change in temperature within 10 % per step.

Despite the unrealistically large extent of the cooling region in some cases,
the hydrogen-line fluxes should be relatively accurate within the other model
assumptions. This is because the hydrogen-line emission occurs in a spatially
thin region (as in the example in Figure 1). Where exactly this region is
located (i.e., possibly at too large a depth) is inconsequential for its
emission properties. The presence of the heated photosphere (as in Calvet &
Gullbring 1998) is more likely to affect the thermal structure of the
postshock region as well as the continuum emission but exploring this is
beyond the scope of this work.

###  6.5 Explaining the non-detections

Recent searches (Cugno et al., 2019; Zurlo et al., 2020; Xie et al., 2020) for
accreting planets in H $\alpha$ have returned non-detections, with 5-$\sigma$
sensitivities down to line-integrated luminosities $L_{\mathrm{H}\,\alpha}\sim
10^{-7}$–$10^{-6}~{}{L_{\odot}}$ beyond $\sim 15$ au ($\gtrsim 100$ mas).
Using the Rigliaco et al. (2012) correlation between $L_{\mathrm{H}\,\alpha}$
and $L_{\mathrm{acc}}$ derived for stars and fixing
$M_{\mathrm{p}}=5~{}M_{\mathrm{J}}$ and $R_{\mathrm{p}}=1.5~{}R_{\mathrm{J}}$,
Zurlo et al. (2020) converted the $L_{\mathrm{H}\,\alpha}$ non-detections to
accretion rate upper limits $\dot{M}\lesssim
10^{-10}$–$10^{-9}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ for most systems. For
those luminosity upper limits, our model for
$M_{\mathrm{p}}=5~{}M_{\mathrm{J}}$ implies instead $\dot{M}\lesssim
10^{-8}$–$10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ (see Figure 4). This is
thus less constraining by a factor of $\sim 100$, which is due to the lower
conversion efficiency from accretion energy to H $\alpha$ in the planetary
case compared to YSOs as commented by Zurlo et al. (2020, see below and ).

The new upper limits of $\dot{M}\lesssim
10^{-8}$–$10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ are still much lower
than the average mass accretion rate of $\dot{M}\sim
10^{-4.5}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ that accumulates a Jupiter mass
within a reasonable disk lifetime (Haisch et al., 2001; see brief review in
Silverberg et al. 2020). A time-dependent accretion rate can solve this
apparent tension. Namely, giant planets could almost reach their final mass
with a much short timescale than the disk lifetime, which means their
observational probability is much lower than that of the disk detection. This
is a feature of runaway gas accretion (e.g., Pollack et al., 1996; Mordasini
et al., 2017), and holds also for gravitational instability (e.g., Greaves &
Rice, 2010; Küffmeier et al., 2017; Nixon et al., 2018; Manara et al., 2018;
Concha-Ramírez et al., 2021; Alves et al., 2020; Segura-Cox et al., 2020;
Schib et al., 2021). Later on, $\dot{M}$ could remain low (Tanigawa & Ikoma,
2007; Tanigawa & Tanaka, 2016). Also, alternating periods of high and low
$\dot{M}$ could occur during the main phase (episodic accretion; Lubow &
Martin, 2012; Brittain et al., 2020; Martin et al., 2021).

Another possibility is that the accretion rate $\dot{M}$ is high but only a
small fraction of at most 1–10 % (the inverse of the 1–2 dex difference in
$L_{\mathrm{H}\,\alpha}$ quoted above) undergoes an accretion shock with
sufficiently high preshock velocity $v_{0}\gtrsim
30~{}\mathrm{km}\,\mathrm{s^{-1}}$ to generate H $\alpha$ (Section A.2). This
could occur if, as in Section 5.2, most gas hits the (thick) CPD relatively
far from the planet, perhaps due to angular momentum effects, and there is no
surface shock. This was studied in Aoyama et al. (2018). Also, extinction can
make the H $\alpha$ flux at the detector lower than the emitted one. This
leads to an underestimate of $\dot{M}$ in some geometry when the PPD, CPD, or
accreting material interrupts the line of sight (e.g., Marleau et al., subm.).

Finally, one obvious explanation is that there are in fact no companions
outside of the inner working angle of detectors up to now ($\gtrsim 100$–200
mas; Close, 2020; Zurlo et al., 2020) as mentioned by Zurlo et al. (2020).
From direct-imaging surveys, several-$M_{\mathrm{J}}$ objects at tens of au or
more are known to be intrinsically rare (e.g., Bowler, 2016; Nielsen et al.,
2019; Vigan et al., 2020). The absence of distant gas giants is also predicted
by classical core accretion theories (e.g., Ikoma et al., 2000; Thommes et
al., 2008). This would leave the $\dot{M}$ of forming planets unconstrained by
the H $\alpha$ survey results.

##  7 Summary and conclusions

Motivated by recent detections of accretion signatures at young planets or
very-low-mass objects (Keppler et al., 2018; Wagner et al., 2018; Haffert et
al., 2019; Eriksson et al., 2020), we have extended the NLTE shock emission
model of Aoyama et al. (2018) to the case that only the planet surface, as
opposed to the circumplanetary disc, is the origin of the hydrogen lines. By
combining the shock spectrum with models for the photospheric emission, we
predict global SEDs (UV to IR) of accreting planets. A possible contribution
from the CPD, relevant at far IR wavelengths (e.g., Zhu, 2015), is not
included. Extinction by the accreting material or the CPD or PPD is neglected
in this work, which we argued (Section 2.4) is a relevant case on its own and
allows us to deal with the complex issue of extinction separately in Marleau
et al. (subm.).

The formation-relevant input parameters of our model are the accretion rate
$\dot{M}$, the planet mass $M_{\mathrm{p}}$, the planet radius
$R_{\mathrm{p}}$, and the filling factor of the accreting region
$f_{\mathrm{fill}}$. To provide guidance, we have fit the radius of forming
planets as a simple but non-monotonic function of $(\dot{M},M_{\mathrm{p}})$
from the results of detailed planet structure calculations (Section A.1.1).
While the structure calculations are not definitive, the fit is an improvement
over using a constant radius as is often done. Other
$R_{\mathrm{p}}(\dot{M},M_{\mathrm{p}})$ relations could be used.

The photospheric effective temperature ${T_{\textnormal{eff}}}$ was derived
approximately self-consistently777In the “Suite of Tools to Model Observations
of accRetIng planeTZ” (St-Moritz) at https://github.com/gabrielastro/St-
Moritz, the $R_{\mathrm{p}}$ and ${T_{\textnormal{eff}}}$ functions are
implemented. from the energy transport from the shock model (Section A.1.2).
Because the shock heats the planet, ${T_{\textnormal{eff}}}$ can easily reach
3000–5000 K for the chosen ranges of $\dot{M}$ and $M_{\mathrm{p}}$ values.

Fixing $f_{\mathrm{fill}}=1$, we have scanned the large parameter space and
shown global SEDs of forming planets (Figure 3). We have also displayed line
luminosities as a function of $\dot{M}$ and $M_{\mathrm{p}}$ for individual
hydrogen lines, focusing on H $\alpha$ as well as H $\beta$, Pa $\alpha$, Pa
$\beta$, Pa $\gamma$, Br $\alpha$, and Br $\gamma$ (Figures 4 and 7) and
discussing their ratios (Figure 8). The data for these and other hydrogen
lines are available upon request.

Our main findings are the following:

  1. 1.

Despite the high ${T_{\textnormal{eff}}}$ of the planet heated by the shock,
the shock contribution to the narrow and broad H $\alpha$ filters of SPHERE
and MagAO dominates over the photospheric contribution (Figure 2).

  2. 2.

At the surface of the planet, the Lyman and Balmer series clearly emerge above
the photosphere over all parameter space, and the Paschen continuum is visible
at low accretion rates. Details depend on the ${T_{\textnormal{eff}}}$ fit but
many lines in the Paschen, Brackett, or other series are visible above the hot
photosphere (Figure 3).

  3. 3.

The H $\alpha$ line luminosity as a function of $\dot{M}$ and $M_{\mathrm{p}}$
is a monotonic function of both and makes it possible to constrain mostly
$\dot{M}$ given an observed value (Figure 4). This is one of the key results.
Applying this tool to current detections yields reasonable constraints
(Section 4). For example, the mass accretion rate of PDS 70 b is estimated as
$\dot{M}=(1.1\pm 0.3)\times 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ from
$L_{\mathrm{H}\,\alpha}=(3.3\pm 0.1)\times 10^{-7}{L_{\odot}}$ (Haffert et
al., 2019). The degeneracy between $\dot{M}$ and $M_{\mathrm{p}}$ can be
lifted at sufficiently high resolution, higher than afforded by MUSE (Aoyama &
Ikoma, 2019; Thanathibodee et al., 2019).

  4. 4.

If there is an accretion shock both on the planet surface and on a
circumplanetary disk, the signal is likely to be dominated by the surface-
shock contribution. The two shocks are expected to be spectrally
distinguishable, with the CPD shock narrower. For PDS 70 b, the resolution of
MUSE is not sufficient, but the minimum requirement depends on how sensitive
future instruments are to the wings of the H $\alpha$ line (Figure 6).

  5. 5.

The line luminosity of other transitions such as H $\beta$, Pa $\alpha$, Pa
$\alpha$, Pa $\beta$, Pa $\gamma$, Br $\alpha$, or Br $\gamma$ is a monotonic
function of $\dot{M}$ and $M_{\mathrm{p}}$ (Figure 7). We compare this to
upper limits for PDS 70 b and the TW Hydra disk. The intensity ratios of these
lines to H $\alpha$ range between $10^{-3}$ and $\approx 1$ (Figure 8), and
their measurement can yield some constraints on the amount of extinction
(Hashimoto et al., 2020).

  6. 6.

Recent H $\alpha$ surveys have resulted in non-detections of planets outside
of $\approx 15$ au ($\gtrsim 100$ mas). Applying the extrapolated YSO
$L_{\mathrm{acc}}$–$L_{\mathrm{H}\,\alpha}$ relationship of Rigliaco et al.
(2012) yields upper limits on the instantaneous accretion rate of
$\dot{M}\lesssim 10^{-10}$–$10^{-9}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ (Zurlo
et al., 2020), but using our models made explicitly for planets implies
instead $\dot{M}\lesssim
10^{-8}$–$10^{-7}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$, which is a less strict
constraint, assuming that planets indeed are present are the surveyed stars
(Section 6.5).

We point out that a determination of $\dot{M}$ from accretion tracers yields a
lower limit on the total accretion rate; the rest of the accreting mass could
be joining the planet without a shock at all (e.g., by boundary-layer
accretion) or with a shock with too-low velocity $v_{0}\gtrsim
30~{}\mathrm{km}\,\mathrm{s^{-1}}$ (Sections A.2.1 and 5.2).

In companion papers, we discuss the use of spectrally-resolved line profiles
for inferring the physical parameters of planets (Aoyama & Ikoma, 2019), study
the correlation between $L_{\mathrm{H}\,\alpha}$ and the accretion luminosity
$L_{\mathrm{acc}}$ in the planetary case (Aoyama et al., submitted), and
assess the absorption of the H $\alpha$ flux by the infalling gas and dust
(Marleau et al., subm.). These models are applied to PDS 70 b also in
Hashimoto et al. (2020) and to Delorme 1 (AB) b in Eriksson et al. (2020).

In Aoyama et al. (submitted), we compare our work to other recent models for
the H $\alpha$ emission associated with accreting planets (Thanathibodee et
al., 2019; Szulágyi & Ercolano, 2020). We emphasise that the postshock
electron populations cannot reach equilibrium as the gas cools (Aoyama et al.,
2018). This and other radiative properties invalidate the key assumptions
behind Storey & Hummer (1995), which was developed for a physically very
different context (e.g., planetary nebulae or H ii regions). Thus Storey &
Hummer (1995) is not appropriate for hydrogen-line luminosity predictions for
the planetary accretion shock (Appendix B). For the stellar context too,
Storey & Hummer (1995) has been suggested to not be appropriate (e.g., Edwards
et al., 2013; Rigliaco et al., 2015; Antoniucci et al., 2017), despite its use
in earlier analyses.

Continued searches with existing and upcoming or proposed instruments are
expected not only to reveal more sources, but will hopefully also increase the
number of detected lines. The first group includes the ZIMPOL subsystem of
VLT/SPHERE (Beuzit et al., 2008; Schmid et al., 2018), as well as VLT/MUSE
(Bacon et al., 2010), LBT/MagAO (Close et al., 2014a, b), SCExAO/VAMPIRES
(Uyama et al., 2020), while to the planned instruments belong MagAO-X (Males
et al., 2018; Close et al., 2018; Close, 2020), KPIC (Jovanovic et al., 2019;
Morris et al., 2020), VIS-X (for H $\alpha$ with $R=15,000$; PI: S. Haffert,
priv. comm.), RISTRETTO888See https://zenodo.org/record/3356296. on the VLT
and later possibly on the ELT (with $R>130,000$, possibly up to $R=150,000$,
covering H $\alpha$; PI: Ch. Lovis; Chazelas et al., 2020), NIRSpec999See
https://jwst-docs.stsci.edu/near-infrared-spectrograph. on the James Webb
Space Telescope (JWST; but note the moderate resolution $R\approx 1000$–2700
and likely high demand for time), HARMONI101010See https://harmoni-
elt.physics.ox.ac.uk. on the ELT (a first-light IFU that will cover H $\alpha$
at $R=3000$ and $\lambda=0.8$–2 $\upmu$m at $R=7000$ or $R=17,000$), and
HIRES/ELT (a second-generation spectrometer planned to cover 1–1.8 $\upmu$m at
$R=100,000$; Marconi et al., 2018; Tozzi et al., 2018, E. Oliva 2020, priv.
comm.). Combined with simulations of forming planets, these rich data sets are
poised to help constrain observationally the complex accretion geometry and
ultimately the origin of gas giants.

We wish to pay tribute to France Allard, who passed away unexpectedly in
October 2020. Her world-leading atmospheric models are used widely in the
observational and theoretical communities, and are an important input in this
work too. Her kind nature and expertise will be missed by many. We thank R.
van Boekel, M. Keppler, A. Müller, J. Bouwman, B. Husemann, M. Samland, D.
Homeier, H. M. Schmid, J. Milli, J. Girard, S. Quanz, G. Cugno, T. Stolker, S.
Edwards, L. Venuti, B. Stelzberg, Ch. Rab, C. Manara, N. Turner, W. Béthune,
M. Bonnefoy, S. Kraus, and O. Ernesto for useful discussions and explanations,
insightful questions, and helpful sharing of data. YA was supported by the
Leading Graduate Course for Frontiers of Mathematical Sciences and Physics.
G-DM acknowledges the support of the German Science Foundation (DFG) priority
program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1
and MA 9185/1-1). G-DM and CM acknowledge support from the Swiss National
Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. Parts of this
work have been carried out within the framework of JSPS Core-to-Core Program
“International Network of Planetary Science (Planet2)” and the NCCR PlanetS
supported by the Swiss National Science Foundation. This research has made use
of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/)
supported by the Spanish MINECO through grant AYA2017-84089 (Rodrigo et al.,
2012; Rodrigo & Solano, 2020). This research was supported in part by the
Japanese–German visitor program of the NCCR PlanetS and by the visitor program
of the German Science Foundation (DFG) priority program SPP 1992.

##  Appendix A Detailed model description

###  A.1 Fitting of planetary properties

####  A.1.1 Radius fit

In principle, the main parameter space for our calculations is
$(\dot{M},M_{\mathrm{p}},R_{\mathrm{p}},f_{\mathrm{fill}})$ along with the
choice of “cold-start” or “hot-start” accretion (high or low radiation
efficiency of the accretion energy; Marleau et al., 2017, 2019). Here, we take
for the sake of definiteness the $R_{\mathrm{p}}(\dot{M},M_{\mathrm{p}})$
relations in the cold- and hot-start populations of the Bern model (Alibert et
al., 2005; Mordasini et al., 2012b, a, 2015, 2017), using data from all time
snapshots111111 The data can be visualized at, and downloaded from, the
“Evolution” section of the Data Analysis Centre for Exoplanets (DACE) platform
under https://dace.unige.ch.. We use the populations CD752 (hot) and CD753
(cold), described and analysed in Mordasini et al. (2012b, a, 2017)
(Generation Ib). Recently, the first results from the Generation III
population syntheses of the Bern model were released (Emsenhuber et al.,
2020a, b; Schlecker et al., 2020), which all assume warm accretion. We
verified that the distribution of points in
$(\dot{M},M_{\mathrm{p}},R_{\mathrm{p}})$ space is very similar between the 1-
and the 100-embryo-per-disk simulations NG73 and NG76, respectively, on the
one hand, and CD752 on the other. Two small differences are that the accretion
rates reached are not quite as high as in Generation Ib (note however that
there are fewer synthetic planets in the region of interest), and that in NG76
the radii can be higher at a given $\dot{M}$ and $M_{\mathrm{p}}$, likely due
to interactions between the embryos. Since the radii are overall similar, we
will keep using the Generation Ib populations in order to cover also high
accretion rates, as could be relevant for instance to accretion outbursts
(e.g., Lubow & Martin, 2012; Brittain et al., 2020; Martin et al., 2021).

These planet structure models were calculated assuming that the planet is at
all times convective. Recent work suggests that forming planets may be in fact
in part radiative (Berardo et al., 2017; Berardo & Cumming, 2017; Cumming et
al., 2018) and thus have a different radius. Nevertheless, the relations
$R_{\mathrm{p}}(\dot{M},M_{\mathrm{p}})$ from the population syntheses provide
a reasonable bracket and reduce the dimensionality of the large parameter
space $(\dot{M},M_{\mathrm{p}},R_{\mathrm{p}},\eta,f_{\mathrm{fill}},\ldots)$.
However, it is clear that this is not meant as a final answer and one could
repeat this study with for example the
$R_{\mathrm{p}}(\dot{M},M_{\mathrm{p}})$ relations of Ginzburg & Chiang
(2019), who find (much) larger radii at a given mass.

After some experimentation, we arrived at the following relatively simple form
for the fitting function121212For convenience, it is provided in different
languages along with the fit coefficients in the “Suite of Tools to Model
Observations of accRetIng planeTZ” (St-Moritz) at
https://github.com/gabrielastro/St-Moritz.:

$\displaystyle R_{1}(\dot{M},M_{\mathrm{p}})=$
$\displaystyle~{}a_{0}+b_{0}{\mathrm{lg}\dot{M}_{2}}+c_{0}\mathrm{e}^{d_{0}{\mathrm{lg}\dot{M}_{2}}}$
$\displaystyle+\left(a_{1}+b_{1}{\mathrm{lg}\dot{M}_{2}}+c_{1}\mathrm{e}^{d_{1}{\mathrm{lg}\dot{M}_{2}}}\right)\left(M_{1}-1\right)$
$\displaystyle+\left(a_{2}+b_{2}{\mathrm{lg}\dot{M}_{2}}+c_{2}\mathrm{e}^{d_{2}{\mathrm{lg}\dot{M}_{2}}}\right)\left(M_{1}-1\right)^{2},$
(A1)

where $R_{1}\equiv R_{\mathrm{p}}/R_{\mathrm{J}}$, $M_{1}\equiv
M_{\mathrm{p}}/M_{\mathrm{J}}$, and
${\mathrm{lg}\dot{M}_{2}}\equiv\log_{10}(\dot{M}/10^{-2}~{}M_{\mathrm{E}}\,\mathrm{yr}^{-1})$.
The fits were performed through gnuplot’s built-in fit routine. Only planets
with $\dot{M}>10^{-5}~{}M_{\mathrm{E}}\,\mathrm{yr}^{-1}$,
$1~{}M_{\mathrm{J}}<M_{\mathrm{p}}<20~{}M_{\mathrm{J}}$, and
$R_{\mathrm{p}}<4~{}R_{\mathrm{J}}$ were used to obtain the fits. We used for
each planet a statistical weight inversely proportional to its radius to have
a more accurate fit at lower radii, for which $v_{0}$ is higher and thus the
accretion signatures a priori stronger. The coefficients for the cold-nominal
population are

⬇

a0 = 1.53; b0 = 0.111;

c0 = 1.06; d0 = 0.906;

a1 = -0.195; b1 = -0.0307;

c1 = 0.0977; d1 = 0.000695;

a2 = -0.250; b2 = 0.000276;

c2 = 0.254; d2 = 0.000214

and for the warm population the coefficients are

⬇

a0 = 0.411; b0 = -0.244;

c0 = 3.45; d0 = 0.762;

a1 = -0.489; b1 = -0.0961;

c1 = 0.652; d1 = 0.353;

a2 = -0.228; b2 = -0.00106;

c2 = 0.226; d2 = 0.000220.

We verified that excluding from the fitting procedure the planets for which
the accretion radius ${R_{\textnormal{acc}}}<10R_{\mathrm{p}}$ (this concerns
only a small fraction of the planets) changed neither the relationships nor
the quality of the fit significantly. The accretion radius
${R_{\textnormal{acc}}}$ is a spherically-averaged estimate of the typical
distance from which the gas is effectively falling onto the planet. It is
defined through

$\frac{1}{{R_{\textnormal{acc}}}}=\frac{1}{{R_{\textnormal{Bondi}}}}+\frac{1}{{k_{\textnormal{Liss}}}{R_{\textnormal{Hill}}}},$
(A2)

where ${R_{\textnormal{Bondi}}}$ and ${R_{\textnormal{Hill}}}$ are the Bondi
and Hill radii and ${k_{\textnormal{Liss}}}=1/4$ (and not, as written in
Mordasini et al., 2012b, ${k_{\textnormal{Liss}}}=1/3$).

The resulting relations $R_{\mathrm{p}}(\dot{M},M_{\mathrm{p}})$ for the cold-
nominal and the warm populations are shown in Figure 9. At high accretion
rates, $R_{\mathrm{p}}$ reaches $\approx 5~{}R_{\mathrm{J}}$ ($\approx
3~{}R_{\mathrm{J}}$) in the “warm” (“cold-nominal”) population. As a function
of mass, the radius monotonically increasing with mass for the “warm
population” and has a minimum near 10–15 $M_{\mathrm{J}}$ for the “cold-
nominal” population.

The fit is overall excellent, with a match to roughly 10 %. For the cold-
nominal population, which displays the largest deviations: Only for masses
$M_{\mathrm{p}}\gtrsim 15~{}M_{\mathrm{J}}$ near
$\dot{M}=10^{-4}$–$10^{-3}~{}M_{\mathrm{E}}\,\mathrm{yr}^{-1}$ is the function
too small, by at most only $\approx 30\,\%$, and at $\dot{M}\gtrsim 3\times
10^{-2}~{}M_{\mathrm{E}}\,\mathrm{yr}^{-1}$ for $M_{\mathrm{p}}$ between 1 and
10 $M_{\mathrm{J}}$ larger or smaller by at most $\approx 30\,\%$. This
reflects in part the intrinsic scatter in the population synthesis results.
For the warm population, the fitted function yields radii also at most 30 %
too small but only towards high masses and low accretion rates. At lower
accretion rates and, in the cold-nominal population, for lower masses than
shown, the fit re-increases but this matches rather well overall the data (not
shown).

Figure 9:  Dependence of the planet radius $R_{\mathrm{p}}$ on accretion rate
$\dot{M}$ and planet mass $M_{\mathrm{p}}$ for planets forming by core
accretion, for the warm- (cold-)population on the left (right). The population
syntheses of Mordasini et al. (2012b) are used but only planets with
${R_{\textnormal{acc}}}>20R_{\mathrm{p}}$ are shown. In the top (bottom) row,
six masses (accretion rates) are considered (see respective legends). The
solid lines display the simple but accurate function of Equation (A.1.1) with
the fitted coefficient as given in the text. Note the logarithmic vertical
scale.

####  A.1.2 Effective temperature

For the photospheric temperature ${T_{\textnormal{eff}}}$, we adopt a semi-
analytical prescription that ensures that, approximately, the total outgoing
flux (shock and photosphere) is equal to the sum of the internal and incoming
energy flux.

At the shock, a portion of the incoming energy is converted into hydrogen-line
and recombination-continua emission. The remaining portion travels downward
into the atmosphere, where it is expected to be thermalised because most of
the energy is in Ly $\alpha$, which can easily be thermalised. What matters
for the structure of the planet is whether this radiation goes deep into the
planet, thereby heating it up, or whether only the top layers are heated up
and re-emit the radiation, on a timescale that is short compared to the
cooling time of the planet. The former outcome corresponds to the “Hot start
(accreting)” case of Mordasini et al. (2012b, their Equation (13)), and the
latter to their “Cold start” case. In both cases, the total luminosity just
outside the planet131313This does not consider the energy recycling in the
accretion flow discussed by Marleau et al. (2017, 2019), and ignores any
absorption in the layers closest to the planet; we take $\eta^{\rm kin}=1$ for
this discussion. is $L\approx L_{\mathrm{int}}+L_{\mathrm{acc}}$, where
$L_{\mathrm{int}}$ is the energy coming from the deep interior.

However, the spectrum of the accreting region (covering a fraction
$f_{\mathrm{fill}}$ of the surface of the planet) differs between the two
extreme cases, as noted by Mordasini et al. (2012b). In the “Hot start”
extreme case the spectrum is entirely thermalised, as given by an atmospheric
model with
${T_{\textnormal{eff}}}^{4}={T_{\textnormal{int}}}^{4}+T_{\mathrm{acc}}^{4}$,
where the “accretion temperature”

$T_{\mathrm{acc}}=\left(\frac{L_{\mathrm{acc}}}{4\pi
R_{\mathrm{p}}^{2}f_{\mathrm{fill}}\sigma}\right)^{1/4},$ (A3)

with $\sigma$ the Stefan–Boltzmann constant. Then there would be no emission-
line spectrum as we have been computing in this work. At the other extreme, in
the “Cold start” case the spectrum from the accreting region is given by the
sum of an atmosphere at ${T_{\textnormal{eff}}}={T_{\textnormal{int}}}$ (i.e.,
not heated up at all by the shock) and the shock emission. Clearly both are
limiting cases. In either case the remaining fraction $(1-f_{\mathrm{fill}})$
of the planet surface has ${T_{\textnormal{eff}}}={T_{\textnormal{int}}}$, and
the global spectral appearance will be a mixture of the spectrum from the
accreting and the non-accreting components, with the proportion set by the
(rotational-phase-dependent) viewing geometry. We set $f_{\mathrm{fill}}=1$ in
this work for simplicity.

Which scenario is likely more accurate? The single-stream, frequency-averaged
simulations of Marleau et al. (2017, 2019) suggest that on the net, only a
small fraction of the incoming $L_{\mathrm{acc}}$ will go in deeper, but that
this small portion is likely (much) higher than the internal luminosity in the
extreme cold starts by Marley et al. (2007). However, for the two-stream,
frequency-dependent calculations presented here, Aoyama et al. (2018)
mentioned that roughly one half of the radiation goes down and one half goes
up, which holds in the limit that the emitting region is optically thin.
Because the radiation transport is more detailed in the models of Aoyama et
al. (2018), we will follow their results. Therefore, we will quantify this
fraction more precisely and use it to derive ${T_{\textnormal{eff}}}$.

Figure 10: Fraction $f_{\mathrm{down}}$ of the incoming kinetic-energy flux
that is in the downward-travelling radiation field at the bottom of the
atmosphere. This includes all hydrogen lines and continua. The solid (dashed)
lines use the hot (cold) population radius fits, and for
$f_{\mathrm{down}}(n_{0},v_{0})$ we use the fit from Equation (A4).

Figure 10 shows the fraction $f_{\mathrm{down}}$ of the incoming kinetic
energy flux that is present in the downward-moving radiation field at the
bottom of the computational domain, where $T$ reaches $10^{4}$ K. For this, we
have first written

$\displaystyle f_{\mathrm{down}}^{\prime}=$
$\displaystyle~{}a_{0}+b_{0}{\mathrm{lg}n_{12}}+c_{0}{\mathrm{lg}n_{12}}^{2}$
$\displaystyle+\left(a_{1}+b_{1}{\mathrm{lg}n_{12}}+c_{1}{\mathrm{lg}n_{12}}^{2}\right)\left(v_{0}-100~{}\mathrm{km}\,\mathrm{s^{-1}}\right)$
$\displaystyle+a_{2}\left(v_{0}-100~{}\mathrm{km}\,\mathrm{s^{-1}}\right)^{2},$
(A4a) $\displaystyle f_{\mathrm{down}}=$
$\displaystyle~{}\min\left(\max\left(f_{\mathrm{down}}^{\prime},0\right),1\right),$
(A4b)

where
${\mathrm{lg}n_{12}}\equiv\log_{10}\left(n_{0}/10^{12}~{}\mathrm{cm}^{-3}\right)$.
Equation (A4a) is a simple biquadratic function in the natural parameters
$v_{0}$ and $\log(n_{0})$ without the $b_{2}$ and $c_{2}$ terms (see below).
Equation (A4b) ensures that $f_{\mathrm{down}}$ remains between 0 and 1 and is
needed only for a small part of the parameter space, for some models on the
grid edge. Using gnuplot’s built-in fit command yielded141414This too is
provided in different languages in the “Suite of Tools to Model Observations
of accRetIng planeTZ” (St-Moritz) at https://github.com/gabrielastro/St-
Moritz.

⬇

a0 = 0.703752; a1 = -0.00527886;

b0 = -0.0967987; b1 = -0.00146833;

c0 = -0.0254579; c1 = -0.000321504;

a2 = -9.91492e-06,

which matches very well the model data (not shown). We had included at first
${\mathrm{lg}n_{12}}$ terms in the $v_{0}^{2}$ term but their coefficients
($b_{2}$ and $c_{2}$) were consistent with zero. We therefore repeated the fit
with the form of Equation (A4a), yielding the reported coefficients. The
result of this fit depends only on our grid of models in $(n_{0},v_{0})$ space
and is thus independent of the population. Next, we used Equations (A7) and
(A8), which relate $(n_{0},v_{0})$ and the macrophysical parameters
$(\dot{M},M_{\mathrm{p}},R_{\mathrm{p}},f_{\mathrm{fill}})$ and are presented
below, along with the radius fits (Equation (A.1.1)) to obtain
$f_{\mathrm{down}}(\dot{M},M_{\mathrm{p}})$.

We find that the fraction $f_{\mathrm{down}}$ does vary at low masses but that
it covers mainly $f_{\mathrm{down}}\approx 0.3$–0.8 between the border of the
brown-dwarf region and a few $M_{\mathrm{J}}$. The fraction depends only
relatively weakly on $\dot{M}$ and much more on $M_{\mathrm{p}}$ (through the
$v_{0}$ dependence), reaching $f_{\mathrm{down}}\approx 1$ at
$M_{\mathrm{p}}\approx M_{\mathrm{J}}$. Note that $f_{\mathrm{down}}$ is
affected to some extent by the non-inclusion of low-temperature ($T\lesssim
10^{4}$ K) cooling processes at the bottom of the computation domain, for
instance from molecules. Currently they are not included, so that in a future
model iteration $f_{\mathrm{down}}$ could be different, once the inclusion of
helium and metals will accelerate the cooling and thus make it computationally
feasible to let the simulations cool down to lower temperatures than
$T=10^{4}$ K. Nevertheless, Figure 10 already provides some guidance.

Having computed $f_{\mathrm{down}}$, we use it to write the photospheric
temperature of the accreting region as

$\displaystyle{T_{\textnormal{eff}}}^{4}$
$\displaystyle={T_{\textnormal{int}}}^{4}+f_{\mathrm{down}}T_{\mathrm{acc}}^{4},$
(A5a) $\displaystyle\Rightarrow L_{\mathrm{phot}}+L_{\mathrm{shock}}$
$\displaystyle=L_{\mathrm{int}}+L_{\mathrm{acc}},$ (A5b)

where $L_{\mathrm{phot}}$ is the luminosity of the photosphere at
${T_{\textnormal{eff}}}$ and $L_{\mathrm{shock}}$ is the total shock
luminosity. In other words, we assume that the component
$F_{\mathrm{down}}=f_{\mathrm{down}}\sigma T_{\mathrm{acc}}^{4}$ from the
shock (lines and recombination continua) is thermalised and re-emitted.
Equation (A5a) ensures that the total upward-travelling radiative flux from
our combined models is
$F=\sigma{T_{\textnormal{eff}}}^{4}+(1-f_{\mathrm{down}})\sigma
T_{\mathrm{acc}}^{4}=\sigma{T_{\textnormal{int}}}^{4}+\sigma
T_{\mathrm{acc}}^{4}$, with the shock flux
$F_{\mathrm{shock}}=(1-f_{\mathrm{down}})\sigma T_{\mathrm{acc}}^{4}$
contained in the models that are the main subject of this work (see Equation
(A5b)). In Equation (A4), to avoid ${T_{\textnormal{eff}}}=0$ K when
$f_{\mathrm{down}}=0$, which occurs at large $(\dot{M},M_{\mathrm{p}})$ in the
cold-start case, we set somewhat arbitrarily ${T_{\textnormal{int}}}=1000$ K.
We verified that this barely affects the function
${T_{\textnormal{eff}}}(\dot{M},M_{\mathrm{p}})$. There is a slight effect at
$\dot{M}\lesssim 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ but this region
of parameter space is of lesser interest.

Figure 11 shows the resulting ${T_{\textnormal{eff}}}$ from Equations (A4) and
(A5a) for both populations for $f_{\mathrm{fill}}=1$. For the warm-population
fit, due to the large radii, ${T_{\textnormal{eff}}}$ ranges only up to
$\approx 2500$ K, whereas in the cold-population ${T_{\textnormal{eff}}}$ can
reach up to 4000 K. Since a logarithmic scale is used for $\dot{M}$ but a
linear one for $M_{\mathrm{p}}$, the dependence of ${T_{\textnormal{eff}}}$ on
$\dot{M}$ appears stronger on $\dot{M}$ than on $M_{\mathrm{p}}$. The low
${T_{\textnormal{eff}}}\approx 1000$–2000 K for the cold-population fit at
large mass $M_{\mathrm{p}}\approx 20~{}M_{\mathrm{J}}$ and radius is a
consequence of $f_{\mathrm{down}}$ going to zero there, with automatically a
large sensitivity on the choice of the fit. Thus it should not be taken very
seriously, and is not of major concern anyway since this part of parameter is
of lesser interest to the present work.

Figure 11:  Approximate effective temperature of accreting planets (Equations
(A4) and (A5a), using ${T_{\textnormal{int}}}=1000$ K; black contours). This
is compared to contours of constant $T_{\mathrm{acc}}$ (red dashed contours).
We fix $f_{\mathrm{fill}}=1$.

We comment briefly on an implicit assumption we make. When we use standard
atmospheric models with Equation (A5a), we are assuming that the emission from
the heated photosphere can be described by that of an isolated object, i.e.,
that if the shock changes the pressure–temperature structure of the atmosphere
(as it is likely), the resulting spectrum is not entirely different. This may
be a strong simplification, and it would be interesting to explore with
dedicated radiative transfer calculations how this modifies the spectral
shape.

Spectra for the BT-Settl models we use (CIFIST2011_2015) exist only for
${T_{\textnormal{eff}}}\geqslant 1200$ K. Therefore, for the case
${T_{\textnormal{eff}}}<1200$ K, we use $1200$ K.

###  A.2 Shock model

####  A.2.1 Shock parameter space

Locally, a strong shock converts most of the mechanical energy into thermal
energy, and the gas temperature increases by orders of magnitude compared to
the preshock value. However, since the temperature much exceeds the radiative
equilibrium temperature at that point, the shock-heated gas in the postshock
region cools rapidly (compared to the postshock flow time) by emitting
radiation.

Given the thinness of the radiation-emitting shock layer compared to the
planetary radius, we assume the shock on the planetary surface to be one-
dimensional as in Aoyama et al. (2018). Then, the shock structure is mainly
determined by two input parameters, namely the hydrogen proton number density
$n_{0}$ and the velocity $v_{0}$ before the shock. The preshock temperature of
the gas $T_{0}$ hardly affects the shock properties, as we discuss below.

As in Aoyama et al. (2018) and e.g. Shapiro & Kang (1987) or Kwan & Fischer
(2011), $n_{0}$ is defined as the immediate preshock number density of
hydrogen nucleons (protons) (i.e., contained in H2, H i, and H+ taken
together). Thus, contrary to the definition, common in the stellar-structure
literature (e.g. Hansen et al., 2004), of $n$ as the number density of all
particles, our $n_{0}$ is independent of the dissociation and ionization
degrees of hydrogen. This definition implies that $n_{0}$ is related to the
(total) preshock gas mass density $\rho_{0}$ by

$X\rho_{0}=n_{0}m_{\mathrm{H}},$ (A6)

where $X$ is the hydrogen mass fraction and $m_{\mathrm{H}}$ is the mass of a
hydrogen atom. We use number ratios given by
$\textrm{H}:\textrm{He}:\textrm{C}:\textrm{O}=1:10^{-1.07}:10^{-3.48}:10^{-3.18}$
(Cox, 2000). As in Aoyama et al. (2018) we do not consider the other elements,
which are neglible. Thus the hydrogen, helium, and metal mass fractions are
respectively $X=0.738$, $Y=0.251$, and $Z=1-X-Y=0.011$.

In the chemistry module, the abundances of the species are defined by $y_{i}$
(see also Iida et al., 2001). The quantity $y_{i}$ is the relative abundance
(in number) of species $i$ (of any particle) with respect to the number of
hydrogen protons. For example, pure H2 has $y_{\mathrm{H}_{2}}=0.5$. Defining
the total $y_{\mathrm{t}}=\sum{y_{i}}$, the usual total number density of all
particles is given by $n=y_{\mathrm{t}}n_{0}$; the number density of particles
of species $i$ is $n_{i}=y_{i}n_{0}$. With the mean weight per particle given
by $\mu_{0}=\sum y_{i}m_{i}/y_{\mathrm{t}}$ for particle masses
$m_{i}m_{\mathrm{H}}$, we have that
$\rho_{0}=y_{\mathrm{t}}\mu_{0}m_{\mathrm{H}}n_{0}$, implying
$X=1/(y_{\mathrm{t}}\mu_{0})$.

The shock input parameters $(n_{0},v_{0})$ are related to the macrophysical,
planet formation parameters
$(\dot{M},M_{\mathrm{p}},R_{\mathrm{p}},f_{\mathrm{fill}})$ by

$\displaystyle v_{0}$
$\displaystyle=\sqrt{\frac{2GM_{\mathrm{p}}}{R_{\mathrm{p}}}}$ (A7)
$\displaystyle n_{0}$ $\displaystyle=\frac{X\dot{M}}{4\pi
R_{\mathrm{p}}^{2}f_{\mathrm{fill}}m_{\mathrm{H}}v_{0}}$ (A8)
$\displaystyle=\frac{X\dot{M}}{\sqrt{32G}\pi
m_{\mathrm{H}}f_{\mathrm{fill}}\sqrt{M_{\mathrm{p}}R_{\mathrm{p}}^{3}}},$ (A9)

where $G$ is the gravitational constant, $M_{\mathrm{p}}$ is the planet mass,
$R_{\mathrm{p}}$ is the planet radius, $\dot{M}$ is the accretion rate, and
$f_{\mathrm{fill}}$ is the filling factor of the shock on the planet surface.
Equations (A7)–(A9) are valid in the limit that the accretion radius
${R_{\textnormal{acc}}}\gg R_{\mathrm{p}}$, with the gas free-falling from
${R_{\textnormal{acc}}}\sim{R_{\textnormal{Hill}}}$ (Equation (A2);
Bodenheimer et al., 2000). Especially at low masses $M_{\mathrm{p}}\lesssim
1~{}M_{\mathrm{J}}$, the accreting gas could be falling in at less than the
free-fall velocity, depending on the thermodynamics and the angular momentum
conservation (Béthune, 2019; Schulik et al., 2019, 2020, but note that the
finite smoothing length of the latter two works imply that converged results
have not quite been reached yet). However, the difference should be small, and
for the classical 1D models that we use as an approximation, the limit
${R_{\textnormal{acc}}}\gg R_{\mathrm{p}}$ usually holds for
$M_{\mathrm{p}}\gtrsim 1~{}M_{\mathrm{J}}$, especially for planets forming at
large distances. Similarly, our model is also applicable when the gas falls
from the inner edge of the CPD, but, in such a case, the $v_{0}$ and the
estimated $L_{\mathrm{H}\,\alpha}$ is smaller by a factor of a few than the
results in this paper. Inserting typical values for stars and planets, the
typical preshock number density is larger in the planetary than in the stellar
case by about a factor of 100 (Zhu, 2015, see Eq. (15) in) but one should keep
in mind that the parameter space is large.

The assumption here is that all of the accreting gas is available for a shock,
whether this turns out to produce H $\alpha$ or not. In reality, some fraction
of the accreting gas could be added to the planet through boundary-layer
accretion (BLA; e.g., Kenyon & Hartmann, 1987; Kley, 1989; Dong et al., 2020),
which does not feature supersonic radial velocities. In this scenario, the
temperature in the boundary layer would not be high enough for H $\alpha$ to
be emitted. Thus converting an observed H $\alpha$ luminosity to planetary
parameters such as accretion rate and mass needs to assume something about the
fraction of the incoming gas that is able to produce H $\alpha$. Put
differently, a measured H $\alpha$ luminosity yields an estimate of the H
$\alpha$-emitting accretion rate, while the total accretion rate could be
higher. However, it is not clear how likely BLA is in the planetary case; Owen
& Menou (2016) argue for BLA but this is based, through the Christensen et al.
(2009) scaling, on magnetic field strengths appropriate of old and faint
planets ($B\approx 0.03$–0.06 kG), not forming nor young, high-luminosity
(Mordasini et al., 2017) objects, which could have $B\sim 1$ kG (Katarzyński
et al., 2016). Thus, assuming that all the accreting gas can undergo a shock
seems reasonable, but high-resolution studies are required to help settle the
question.

Figure 12:  Preshock velocities of forming planets. Left panel: Preshock
velocity against planet radius or radial distance (for a shock on the CPD).
Bern cold-nominal population synthesis planets with
${R_{\textnormal{acc}}}\geqslant 10R_{\mathrm{p}}$ are shown, coloured by
$\dot{M}$, and the free-fall velocity (Equation (A7)) is drawn for different
planet masses (dashed lines). Pink crosses are from the Alcalá et al. (2017)
YSO sample. Below $v_{0}\approx 30~{}\mathrm{km}\,\mathrm{s^{-1}}$, there is
no line emission (H $\alpha$, H $\beta$, etc.) if the gas comes in molecular
(gray region). Right panel: Velocity against mass density or hydrogen proton
number density (lower, upper axes; Equation (A6)). Population synthesis points
are shown, colored by the preshock $\mu_{0}$ (for $\mu_{0}\approx 1.3$, only
half the points are shown). Dashed lines: constant $\log\dot{M}=-7.5$ to
$-4.5$ ($M_{\mathrm{J}}\,\mathrm{yr}^{-1}$) and $R_{\mathrm{p}}=1.5$ (black)
or $3~{}R_{\mathrm{J}}$ (gray). We fix $f_{\mathrm{fill}}=1$.

In Figure 12 we plot the preshock velocity as a function of $\dot{M}$
(colour), mass (to be read off from the dashed lines), and radius planet ($x$
axis) for the cold-nominal Bern population synthesis (Mordasini et al., 2012a,
2017). The hot population objects (not shown) have larger radii but an
otherwise similar distribution. All snapshots are included but we restrict
$\dot{M}\geqslant 10^{-8}~{}M_{\mathrm{J}}\,\mathrm{yr}^{-1}$ and
$M_{\mathrm{p}}\geqslant 0.1~{}M_{\mathrm{J}}$. Planets with
$R_{\mathrm{p}}<10{R_{\textnormal{acc}}}$ are excluded, but this affects only
a small fraction of the points; for example, planets with
$M_{\mathrm{p}}\approx 0.1~{}MJ$ and with a small Hill sphere might be in the
attached phase even at high $\dot{M}$. Essentially all points above 0.5
$M_{\mathrm{J}}$ have $R_{\mathrm{p}}>10{R_{\textnormal{acc}}}$. Thus, $v_{0}$
is indeed given by Equation (A7) to 5 % or better.

Figure 12a shows that typical preshock velocities are $v_{0}\approx 50$–200 km
s-1, and that there is a maximal $v_{0}\approx
170~{}\mathrm{km}\,\mathrm{s^{-1}}$ for masses between $M_{\mathrm{p}}\approx
13$ and $\approx 35$ $M_{\mathrm{J}}$. (For the hot population, it is also
$v_{0}\approx 170~{}\mathrm{km}\,\mathrm{s^{-1}}$ but with a thicker spread of
the horizontal portion there, down to $v_{0}\approx
150~{}\mathrm{km}\,\mathrm{s^{-1}}$.) The points along this locus of
approximately constant $v_{0}$ (which implies that $M_{\mathrm{p}}\propto
R_{\mathrm{p}}$ approximately) are objects which have not yet burned their
deuterium and are all younger than 10 Myr. These points span a range of ages
and masses. Thus their distribution is different from the D-burning “shoulder”
in plots of the radius of an object of a given mass as a function of time, in
which the D burning stalls the cooling (e.g., Mollière & Mordasini, 2012, and
references therein). Also, these points with $M_{\mathrm{p}}\propto
R_{\mathrm{p}}$ are not in conflict with the classical result $M\propto
R^{-1/3}$ in that mass regime (Zapolsky & Salpeter, 1969) since that holds for
“zero-temperature” (i.e., degenerate) objects, whereas the objects here are
still forming and thus hotter, also in the sense of not being completely
degenerate. Older, post-deuterium-burning objects have a smaller radius and
thus a higher $v_{0}$ but they are not plotted because of the selection on
$\dot{M}$. This group of constant-$v_{0}$ points could show up as a pile-up in
a histogram of $v_{0}$ for forming planets.