[orcid=0000-0002-0140-4475] [1] [1] 1]organization=Department of Astronomy,
University of Michigan, city=Ann Arbor, postcode=48109, state=MI, country=USA

[orcid=0000-0002-8167-1767] 2]organization=Department of Physics, University
of Michigan, city=Ann Arbor, postcode=48109, state=MI, country=USA

[cor1]Corresponding author [fn1]Fannie and John Hertz Foundation Fellow

# Formation and Structure of Circumplanetary Disks and Envelopes during the
Late Stages of Giant Planet Formation

Aster G. Taylor<EMAIL_ADDRESS>[  Fred C. Adams  [

###### Abstract

Giant planets are expected to form within circumstellar disks, which shape
their formation history and the local environment. Here, we consider the
formation and structure of circumplanetary disks that arise during the late
stages of giant planet formation. During this phase, when most of the final
mass is accumulated, incoming material enters the Hill sphere and falls toward
the planet. In the absence of torques, the falling parcels of gas conserve
their specific angular momentum and collect into a circumplanetary disk.
Generalizing previous work, we consider a range of possible geometries for the
flow entering the sphere of influence of the planet. Specifically, we consider
five geometric patterns for the inward flow, ranging from concentration toward
the rotational poles of the system to isotropic flow to concentration along
the equatorial plane. For each case, we derive analytic descriptions for the
density field of the infall region, the disk surface density in the absence of
viscosity, and steady-state solutions for viscous disks. These results, in
turn, specify the luminosity contributions of the planet, the circumplanetary
disk, and the envelope. These power sources, in conjunction with the
surrounding material, collectively determine the observational appearance of
the forming planet. We conclude with an approximate determination of these
radiative signatures.

###### keywords:

Planet formation (1241) Protoplanetary disks (1300) Planetary system formation
(1257) Solar system formation (1530) Extrasolar gas giant planets (509)

##  1 Introduction

Most giant planets are thought to form through the core accretion paradigm,
and most of the planetary mass is accreted during the final stage of the
process. Due to conservation of angular momentum, the majority of the incoming
material cannot fall directly onto the planetary surface, but instead collects
into a disk structure surrounding the planet. The formation, structure, and
evolution of these disks thus represents a crucial part of the
planet–formation paradigm. Our goal is to construct an analytic description of
this process by generalizing previous work to include a wider range of
boundary conditions.

The core accretion paradigm can be broken up into three distinct stages
(Pollack et al., 1996). First, rocky and icy material collects into a core
with mass $M_{p}\simeq 10\,M_{\oplus}$. After crossing this threshold, the
protoplanet enters the second stage and gathers a hydrostatically-supported
gaseous envelope extending from the core surface out to (approximately) the
Hill radius. The gaseous envelope will slowly cool and collapse, lengthening
the formation time while adding little mass. After the envelope mass becomes
comparable to that of the rocky core, hydrostatic balance is insufficient to
support the envelope. As a result, pressure no longer supports the infalling
material, and the growing planet begins to gain mass in a runaway process.
Although the second intermediate phase takes the longest, most of the mass of
the planet is gathered during the third and final phase. As incoming material
transitions from the background circumstellar disk into the sphere of
influence of the planet, it retains its specific angular momentum inherited
from the disk (in the absence of strong torques from magnetic fields and other
mechanisms). As a result, the incoming material cannot fall directly to the
planetary surface, but instead collects into a circumplanetary disk. The
structure and evolution of this planet–disk system determine the properties of
the forming planet during the third stage of formation. We determine the
properties of infall regions and circumplanetary disks that arise during this
late stage of giant planet formation.

A great deal of numerical work has been carried out concerning the problem of
giant planet formation in circumstellar disks, including the third stage when
rapid mass accretion occurs (e.g., Szulágyi et al. 2016; Lambrechts et al.
2019, and many others). These simulations often show that circumplanetary
disks form alongside planets, as expected due to conservation of angular
momentum (see also Machida et al. 2008; Szulágyi and Mordasini 2017; Fung et
al. 2019). However, even state-of-the-art numerical treatments remain somewhat
limited in resolution, both temporally (with integration times much shorter
than the typical $\sim 10^{5}$ orbits of the third stage) and spatially (often
using a relatively modest number of grid cells to span the Hill sphere of the
forming planet), thus motivating a complementary analytic treatment. More
significantly, the results of different numerical treatments are not
consistent with each other. Specifically, we are interested in the angular
distribution of material entering the Hill sphere. A set of simulations finds
that “the vast majority of the mass flows into the Hill sphere near the
equator” (Ayliffe and Bate, 2009), while a competing study finds that
“radiative hydrodynamical simulations reveal a steady-state gas flow, which
enters through the poles” (Lambrechts and Lega, 2017). Although these authors
find similar meridional flow patterns, the amount of mass loaded onto a given
streamline shows significant differences in these studies. These differences
include geometrical effects and go beyond the differing outcomes resulting
from the chosen equations of state, which lead to some numerical simulations
finding that disks readily form, whereas others do not. In any case, the
problem of circumplanetary disk formation is not yet fully understood, and a
range geometries for the incoming flow remain possible.

As a result, we here consider how the structure of a forming gas giant and its
circumplanetary disk depend on the geometry of the inflowing material. This
work expands upon an earlier treatment (Adams and Batygin, 2022) that
calculated the theoretical structure of the disk and envelope for the case
where flow enters the Hill sphere isotropically.

In recent years, observational capabilities have progressed to the point that
investigations of the regions immediately surrounding forming gaseous planets
are feasible. Although current observations do not have the angular resolution
to directly image the size scales of the planets themselves (with expected
radii $R_{p}\sim 10^{10}$ cm), they can begin to resolve circumplanetary disks
(with expected sizes of $R_{d}\sim 2\times 10^{12}$ cm $\sim 0.1$ AU),
although such scales are currently at (or near) the resolution limit. The
first tentative detection of a circumplanetary disk around a forming planet
(Benisty et al., 2021) is expected to be the first of many. Additional
observations in the coming years will thus elucidate the formation processes
for giant planets. Corresponding models of circumplanetary structure are
therefore needed. The goal of this work is to provide the first step towards
the requisite theoretical understanding.

This paper is organized as follows. Sec. 2 specifies the environment for
forming planets, including the outer boundary conditions, as determined by the
background circumstellar disk. For each of the inward flow patterns, Sec. 3
finds the corresponding infall solutions, including the density field, the
velocity field, and the column density of the envelope. These quantities
determine the structure and formation of the circumplanetary disks, which are
considered in Sec. 4, with a focus on their surface density distributions. We
present solutions for disks that are built up from the infalling flow with no
further evolution and for the case of steady-state viscous accretion. In Sec.
5, we derive the luminosity of the planet, disk, and envelope, making several
simplifying assumptions, and develop an approximate description of their
spectral energy distributions. Finally, the paper concludes in Sec. 6 with a
summary of our results and a discussion of their implications.

Figure 1: The infall functions used in this paper, which characterize the
asymmetry in the material falling into the Hill sphere. The functions range
from primarily-polar to primarily-equatorial infall. Each function is
normalized such that $\oint f_{i}{\rm d}\Omega=1$.

##  2 The Background Circumstellar Disk Environment

In this paper, we consider the formation of gas giant planets, which actively
accrete material from their parent circumstellar disks and subsequently form
their own circumplanetary disks. The majority of gas giant planets are thought
to form within circumstellar disks hosted by at least one star, and these
disks provide the boundary conditions for the process of planet formation. As
outlined above, after the growing planet reaches a mass threshold $M_{p}\simeq
20$ ${\rm M}_{\oplus}$, it can gain mass rapidly from its nearby environment.

The boundary between the environment of the forming planet and the background
circumstellar disk is nominally provided by the Hill radius
$R_{H}=a(M_{p}/3M_{\ast})^{1/3}$, where $a$ is the semi-major axis of the
planetary orbit. However, if the Bondi radius $R_{B}$ = $GM/c_{s}^{2}$ is
smaller than $R_{H}$ (where $c_{s}$ is the local sound speed), then $R_{B}$
marks the outer boundary of the planetary environment. In addition, if the
scale height of the circumstellar disk $H=c_{s}/\Omega$ is small enough, then
the outer boundary can be modified further by the thinness of the disk. In
practice, the Hill radius, the Bondi radius, and the disk scale height are
roughly comparable, but the Hill radius provides the boundary in most
instances. In this paper, we use the Hill radius as a working approximation
for the location of the boundary, but the results can easily be modified for
other choices. Moreover, we will use the term ‘Hill radius’ as a proxy for
this more general picture of the outer boundary.

In order to determine how material flows toward a nascent planet, and forms a
circumplanetary disk, we must specify the conditions at the outer boundary.
First, we assume that the background circumstellar disk feeds material into
the sphere of influence of the planet with a net rate ${\dot{M}}$. Numerical
simulations show that some material that initially flows into the Hill sphere
will subsequently flow back out. In our formulation of the problem, the mass
accretion rate ${\dot{M}}$ only includes the material that remains within
$R_{H}$. Unless specified otherwise, we use a fiducial (net) mass accretion
rate of 3 $M_{\rm J}$/Myr.111We note that larger accretion rates are sometimes
advocated (e.g, Lubow et al. 1999; D’Angelo et al. 2002; Bate et al. 2003;
Ayliffe and Bate 2009) and be be easily accommodated within the analytic
framework developed herein.

To determine the size of the Hill sphere, we must specify the instantaneous
mass $M_{p}$ of the planet along with its current orbital location $a$. As a
result, the collection of variables $(a,M_{p},{\dot{M}})$ is required to
specify the boundary conditions for the planet formation problem. In addition,
for a protoplanet with given mass $M_{p}$, we must specify its radius $R_{p}$
and the magnetic field strength $B_{p}$ on its surface, which provide inner
boundary conditions. For bodies of roughly Jovian mass, the radius is a slowly
varying function of $M_{p}$, so we take the planetary radius to have a
characteristic value $R_{p}=10^{10}$ cm $\simeq 1.3R_{J}$ (e.g., Marley et al.
2007), although young planets could in principle have radii that are even
larger compared to Jupiter.

When its mass becomes sufficiently large, a forming planet will begin to
affect the background circumstellar disk, potentially creating gaps, spiral
density waves, and additional structure. The planetary mass required for gap
opening depends on the viscosity of the disk and other properties (e.g., see
Fung et al. 2014), but often occurs for masses $M\sim 50-100M_{\oplus}$. When
planets accrete from gaps, the incoming material is no longer azimuthally
symmetric, but rather enters the Hill sphere along specific directions.
Because the gap region is empty, or at least has lower surface density, the
incoming material can enter the Hill sphere along directions more concentrated
towards the equator, as it is accreted from the regions beyond the gap. The
advent of gap opening thus affects both the azimuthal and polar angle
dependence of the incoming material. In this treatment, we retain the
assumption of azimuthal symmetry (although this assumption can be relaxed
using the formulation developed here). In addition to being convenient, this
approximation is justified because the infalling material spirals inward to
the circumplanetary disk, where it rapidly becomes evenly distributed in
azimuth. Note that the orbital timescale within the circumplanetary disk is
$\sim 1$ yr, whereas the infall timescale is $\sim 1$ Myr, so the disk has
ample time to reestablish axial symmetry. On the other hand, the dependence of
the incoming material on the polar angle $\theta$ can play an important role,
as it determines the specific angular momentum and, hence, the location within
the circumplanetary disk where the material lands. As a result, we consider a
range of $\theta$-dependent inflow patterns, as outlined below.

Finally, we must specify the geometry of the material flowing inward through
the outer boundary of the protoplanetary envelope. We consider azimuthal
symmetry and latitudinal asymmetry in the inflow. More specifically, we use
five different inflow functions $f_{i}(\theta_{0})$, which set the mass flux
through the Hill sphere at a polar angle $\theta_{0}$. These functions range
from flow that is concentrated along the poles to isotropic flow to flow that
is concentrated along the equatorial directions. These are written as
trigonometric functions, although the absolute value is always taken, so that
they are mirrored around $\theta_{0}=\pi$ (we will not write the absolute
value signs for convenience, but always assume full symmetry). The flow
functions are normalized so that the total mass inflow is $\dot{M}$, which can
be enforced by setting the integral $\oint f_{i}d\Omega=1$. With the
definition $\mu_{0}\equiv\cos\theta_{0}$, we can write the asymmetry functions
for the inflow in the form

$f_{i}(\mu_{0})=\begin{dcases}3\mu_{0}^{2}&\text{polar},\\\
2\mu_{0}&\text{quasipolar},\\\ 1&\text{isotropic},\\\
\frac{4}{\pi}(1-\mu_{0}^{2})^{1/2}&\text{quasieq.},\\\
\frac{3}{2}(1-\mu_{0}^{2})&\text{equatorial}.\end{dcases}$ (1)

With these specifications, the mass inflow at a point on the Hill sphere is
determined by $\dot{M}f_{i}(\mu_{0})$. The asymmetry functions are shown in
Fig. 1.

Figure 2: Infall trajectories projected into the poloidal plane, extending
from the Hill sphere (dashed line) to the circumplanetary disk (in the
equatorial plane along the $x$-axis). These curves are solutions to the orbit
equation (2). The inflow is roughly spherical at large radii, but becomes more
complex near $x=1\,R_{C}$, $z=0.1\,R_{C}$. The centrifugal radius $R_{C}$ is
the outer edge of the disk.

##  3 Infall-Collapse Solution for the Protoplanetary Envelope

With the boundary conditions and other system properties specified (see Sec.
2), we can determine how incoming material falls inward to toward the planet
and its disk.

###  3.1 The Orbit Equation

We now consider the behavior of matter once it has entered the planetary
sphere of influence. As a leading-order approximation, we neglect pressure
forces and assume that the infalling material travels on zero-energy ballistic
orbits (see, e.g., Ulrich 1976; Chevalier 1983; Cassen and Moosman 1981). As
outlined in Appendix A of Adams and Batygin (2022), under typical conditions,
pressure forces correspond to a $\sim 3-30\%$ correction near the Hill radius
and become increasingly smaller as the material falls inward. In general, the
collapse flow will tend to follow ballistic trajectories as long as cooling is
efficient. In contrast, however, our solutions are not applicable in the
regime of adiabatic collapse, where pressure can dominate over gravitational
forces. In this case, however, the adiabatic envelope does not collapse to
form a planet (Fung et al., 2019), limiting the relevance of this
process.222See also Martin and Lubow (2011) for a comparison of ballistic and
hydrodynamic models of disk inflow.

In the ballistic limit, the gas must conserve its angular momentum while
falling from a point on the sphere given by its initial polar angle
$\theta_{0}$. When the gas parcel enters the Hill sphere, it has an angular
velocity $\Omega$ inherited from the circumstellar disk. Its distance from the
planet–disk rotation axis is given by $R_{H}\sin\theta_{0}$, and so the
specific angular momentum is given by $j=\Omega R_{H}^{2}\sin^{2}\theta_{0}$.
As the parcel moves inwards, it conserves its angular momentum, and thus
follows an orbital equation given by

$1-\frac{\mu}{\mu_{0}}=\zeta(1-\mu_{0}^{2})\,.$ (2)

In Eq. (2), the dimensionless parameter
$\zeta\equiv(\Omega^{2}R_{H}^{4})/(GM_{p}r)$ and $\mu$ = $\cos\theta$, where
$\theta$ is the polar angle in the spherical coordinate system centered on the
planet.333Note that solving this equation directly for $\mu_{0}$ in terms of
($r$, $\mu$) requires the solution to a cubic and is thus complicated for most
of the parameter space, although it is analytically exact. For any given
starting position, the parcel of gas will intersect the disk plane at $\mu=0$,
shock, and join onto the disk. The size of the disk is set by the parcels with
the largest specific angular momentum, which enter the Hill radius from the
equatorial plane and intersect the disk at the centrifugal radius. We can
solve Eq. (2) with $\mu=0$, $\mu_{0}\rightarrow 0$ to find that
$\zeta=R_{C}/r$. Assuming Keplerian orbits around the total mass $M$ that has
fallen at a given epoch, we can also write $\Omega$ in terms of $G$, $M$, and
$R_{H}$ to find that

$R_{C}=\frac{R_{H}}{3}\,.$ (3)

As the planet increases its mass, the Hill radius $R_{H}$ increases
accordingly, and the circumplanetary disk grows so that it always extends to
one third of the Hill radius (Quillen and Trilling, 1998). The trajectories of
the incoming parcels are shown in Fig. 2, where the paths are projected into
the poloidal plane (see also Tanigawa et al. 2012; Adams and Batygin 2022).

Using Eq. (2) and conservation of energy, one can solve for the velocity field
of the collapse flow. After defining a velocity scale

$v_{0}\equiv(GM_{p}/r)^{1/2}\,,$ (4)

we can then write the velocity components in the forms

$\displaystyle v_{r}=$ $\displaystyle-
v_{0}\left[2-\zeta(1-\mu_{0}^{2})\right]^{1/2}$ (5a) $\displaystyle
v_{\theta}=$ $\displaystyle
v_{0}\Big{[}\frac{1-\mu_{0}^{2}}{1-\mu^{2}}(\mu_{0}^{2}-\mu^{2})\zeta\Big{]}^{1/2}$
(5b) $\displaystyle v_{\phi}=$ $\displaystyle
v_{0}\Big{[}\frac{1-\mu_{0}^{2}}{(1-\mu^{2})^{1/2}}\zeta^{1/2}\Big{]}\,.$ (5c)

Note that the variable $\mu_{0}$ is specified by Eq. (2) for any input
location $(r,\mu)$, and so the velocity field is well defined throughout the
envelope. The fact that $v_{r}<0$ indicates that the velocity is radially
inward (see Adams and Batygin 2022; and others).

Table 1: Parameter Values. Fiducial values for model parameters, which are
used as defaults throughout this paper. While variations in these properties
are possible, the analytical nature of our solutions allows the values to be
changed in a straightforward manner.

System Properties
---
$M_{\star}$ | 1 | M⊙
$a$ | 5 | au
$\dot{M}$ | 3 | MJ/Myr

Opacity Properties
---
$\kappa_{0}$ | 10 | cm2 g-1
$\nu_{0}$ | $10^{14}$ | Hz
$\eta$ | 1 |

Planet Properties
---
$M_{p}$ | 1 | MJ
$R_{p}$ | $10^{10}$ | cm
$B_{p}$ | 500 | G

Figure 3: Density distribution of the protoplanetary envelope for each infall
geometry, given by Eq. (7). The colormap and contours are shared by all the
subfigures. The contours range from $10^{-17}$ to $10^{-13}$, separated by
$1/4$ in log space. The parameter values are given by Table 1.

###  3.2 Envelope Density Distribution

The density of material in the envelope is given by the conservation of mass
along streamlines (again see Ulrich 1976 and subsequent references). For any
given point in the envelope $(r,\mu)$, we can solve Eq. (2) to find the
initial point where the parcel crossed the Hill sphere, denoted as
$(R_{H},\mu_{0})$. Conservation of mass then implies that

$\displaystyle\rho(r,\mu)$ $\displaystyle|v_{r}(r,\mu)|\,{\rm d}A(r,\mu)$ (6)
$\displaystyle=$ $\displaystyle\rho(R_{H},\mu_{0})|v_{r}(R_{H},\mu_{0})|\,{\rm
d}A(R_{H},\mu_{0})\,.$

In Eq. (6), ${\rm d}A$ is the area element at a given radius and polar angle.
Since $\rho(R_{H},\mu_{0})=\dot{M}f_{i}(\mu_{0})/|v_{r}(R_{H},\mu)|$ and ${\rm
d}A(r,\mu)=r^{2}{\rm d}\mu{\rm d}\phi$, we find that

$\displaystyle\rho(r,\mu)=$ $\displaystyle\frac{\dot{M}f_{i}(\mu_{0})}{4\pi
r^{2}|v_{r}(r,\mu)|}\frac{{\rm d}\mu_{0}}{{\rm d}\mu}$ (7) $\displaystyle=$
$\displaystyle\frac{\dot{M}f_{i}(\mu_{0})}{4\pi
r^{2}|v_{r}|}[1+\zeta(3\mu_{0}^{2}-1)]^{-1}\,.$

In Eq. (7), we used Eq. (2) to find ${\rm d}\mu_{0}/{\rm d}\mu$. Notice that
this result is independent of the form of the inflow function $f_{i}$, so that
extending this result using Eq. (1) is straightforward.

In Fig. 3, we show the resulting density distributions for the different
inflow functions, with reference parameters given in Table 1. Unless otherwise
stated, these parameter values are used throughout the rest of this paper. The
envelope density distributions show significant differences across the inflow
functions. The majority show an increased density around $r=R_{H}/3$, where
the disk begins and a significant fraction of the infall tracks intersect the
plane (see Fig. 2). However, there are significant differences in the strength
of this and other qualitative features. The isotropic, quasipolar, and polar
inflows show relatively similar density distributions. All three have a second
point of high density near the planet, which is expected as a result of the
relatively high concentration of low-angular-momentum inflow. As the inflow
functions become more concentrated towards the pole, the strength of the near-
planet maximum increases and the strength of the disk-edge maximum decreases,
with the disk-edge maximum nearly disappearing in the fully polar case. In
general, however, these density distributions are significantly more spatially
uniform than in the quasiequatorial and equatorial cases.

For the flow functions that focus material toward equatorial directions,
angular momentum conservation significantly concentrates the density towards
the disk-edge maximum, and relatively little density remains in the direction
of the planetary pole, or even the disk interior. Both of these inflow
distributions then essentially lead to a toroidal circumplanetary envelope,
allowing for optically thin observations of the forming planet and the disk
along the poles. It is also worth noting that these density distributions
require a disk to accrete material onto the planet, since almost no material
falls in via the poles. Although this trend also applies in general, the
effect is not as severe for polar inflow functions.

In this section, we have assumed azimuthal symmetry in the inflow function.
This assumption is justified by the fact that the trajectories into the system
depend primarily on the angular momentum, and thus only on latitude, not
longitude. Due to the difference in timescales, asymmetry in the axial
variable will be readily erased within the structure of both the disk and the
planet. However, asymmetries in the background circumstellar disk, such as
spiral density waves and/or a “headwind” from sub-Keplerian rotation, could
cause the envelope to only partially cover the planet (Cimerman et al., 2017;
Krapp et al., 2022; Kuwahara and Kurokawa, 2024). These effects are expected
to be subdominant, however, and will be important mostly for the density
distribution of the envelope.

###  3.3 Equivalent Spherical Infall

We now find the equivalent spherical infall, which is obtained by averaging
the density distribution over the polar coordinate $\mu$. This average is
given by

$\langle\rho\rangle=\int\displaylimits_{0}^{1}\frac{\dot{M}f_{i}(\mu_{0})}{4\pi
r^{2}|v_{r}|}\frac{{\rm d}\mu_{0}}{{\rm d}\mu}\,{\rm
d}\mu=\int\displaylimits_{\eta}^{1}\frac{\dot{M}f_{i}(\mu_{0})}{4\pi
r^{2}|v_{r}|}\,{\rm d}\mu_{0}\,.$ (8)

The lower limit $\eta$ of the integral depends on the radius $r$, since the
streamlines that strike the disk are irrelevant. For convenience, we can
define a dimensionless radius coordinate $u\equiv r/R_{C}$, and therefore find
that

$\eta=\begin{dcases}\sqrt{1-u}&u<1;\\\ 0&u>1\,.\end{dcases}$ (9)

The lower bound $\eta$ is defined because the streamlines begin to intersect
the disk within $u<1$. As a result, there is a lower limit on the values of
$\mu_{0}$ that contribute to the infall for $u<1$, which can be found by
solving Eq. (2) for points where $\zeta>1$, $\mu=0$.

The density profile that results from Eq. (8) is of the general form

$\langle\rho\rangle=Cr^{-3/2}\mathcal{A}(u)\,.$ (10)

In Eq. (10), the constant $C\equiv\dot{M}_{p}/4\pi\sqrt{2GM_{p}}$ and the
asphericity function $\mathcal{A}(u)$ is defined by

$\mathcal{A}(u)\equiv\sqrt{2}\int\displaylimits_{\eta}^{1}f_{i}(\mu_{0})[2-\zeta(1-\mu_{0}^{2})]^{-1/2}\,{\rm
d}\mu_{0}\,.$ (11)

It is convenient to define

$a^{2}\equiv\frac{2-\zeta}{\zeta}=2u-1\,,$ (12)

so that Eq. (11) becomes

$\mathcal{A}(u)=\sqrt{2u}\int\displaylimits_{\eta}^{1}\frac{f_{i}(\mu_{0})}{[a^{2}+\mu_{0}^{2}]^{1/2}}\,{\rm
d}\mu_{0}\,.$ (13)

We now solve Eq. (13) for the inflow functions given in Eq. (1). For
convenience, and due to repeated patterns in the solutions, we will define a
function

$\mathcal{D}(u)=\begin{dcases}\ln\Big{(}\frac{1+\sqrt{2u}}{\sqrt{1-u}+\sqrt{u}}\Big{)}&u<1\,;\\\
\ln\Big{(}\frac{1+\sqrt{2u}}{\sqrt{2u-1}}\Big{)}&u>1\,.\end{dcases}$ (14)

With this definition, the isotropic case (as given in Adams and Batygin 2022;
see also Adams and Shu 1986), has the solution

$\mathcal{A}(u)=\sqrt{2u}\mathcal{D}(u)\,.$ (15)

For a quasipolar inflow, we find that

$\mathcal{A}(u)=2\sqrt{2u}\begin{dcases}\big{[}\sqrt{2u}-\sqrt{u}\big{]}&u<1\,;\\\
\big{[}\sqrt{2u}-\sqrt{2u-1}\big{]}&u>1.\end{dcases}$ (16)

Polar inflow produces

$\mathcal{A}(u)=\frac{3}{2}\sqrt{2u}\begin{dcases}\begin{aligned}
\big{[}&\sqrt{2u}-(2u-1)\mathcal{D}(u)\\\
&-\sqrt{u(1-u)}\big{]}\end{aligned}&u<1\,;\\\
\big{[}\sqrt{2u}-(2u-1)\mathcal{D}(u)\big{]}&u>1\,.\end{dcases}$ (17)

For the quasiequatorial case, the integral must be expressed in terms of
elliptic integrals. Specifically,

$\mathcal{A}(u)=\frac{4}{\pi}\begin{dcases}\begin{aligned}
2u\big{[}&F\big{(}\sin^{-1}(\sqrt{u}),1/2u\big{)}\\\
&-E\big{(}\sin^{-1}(\sqrt{u}),1/2u\big{)}\big{]}\end{aligned}&u<1\,;\\\
\begin{aligned} \sqrt{2u/(2u-1)}\big{[}2uK(1/(1-2u))\\\
-(2u-1)E(1/(1-2u))\big{]}\end{aligned}&u>1\,.\end{dcases}$ (18)

In Eq. (18), $K(x)$ is the complete elliptic integral of the first kind,
$E(x)$ is the complete elliptic integral of the second kind, $E(x,m)$ is the
incomplete elliptic integral of the second kind, and $F(x,m)$ is the
incomplete elliptic integral of the first kind. Since these can be simply
numerically evaluated, this form is relatively simple to use, although it is
analytically complex.

Figure 4: The equivalent spherical density for each infall function. The
vertical dashed line is at $r=R_{C}$, the centrifugal radius/disk edge. The
spike at this point is a result of a large fraction of material striking the
disk and disappearing from the envelope.

Finally, for the equatorial case, we find that

$\mathcal{A}(u)=\frac{3}{4}\sqrt{2u}\begin{dcases}\begin{aligned}
\big{[}&(2u+1)\mathcal{D}(u)-\sqrt{2u}\\\
&+\sqrt{u(1-u)}\big{]}\end{aligned}&u<1\,;\\\
\big{[}(2u+1)\mathcal{D}(u)-\sqrt{2u}\big{]}&u>1\,.\end{dcases}$ (19)

We show these density profiles in Fig. 4. The behavior of these profiles are
nearly identical beyond the disk, since the same amount of material is falling
through the Hill sphere, and no trajectories have yet impacted the disk. At
$r\leq R_{C}$, the profiles diverge significantly. The equatorial and
quasiequatorial cases exhibit a significant deficit of density, as the
equatorial trajectories impact the disk first (see Fig. 2). The results for
the remaining flow functions show a modest feature near the centrifiugal
radius, along with an shallow increase toward smaller radii. This trend is
most significant for polar inflows, because relatively little mass is accreted
along the equatorial streamlines which strike the disk first.

The density distribution has a steep slope (power-law index $p=3/2$) in the
outer envelope, and a much shallower slope ($p\leq 1/2$) in the inner parts
$(r<R_{C})$. As a result, the contribution to the column density
($N\sim\int\\!\\!\rho\,{\rm d}r$) is concentrated near the centrifugal radius
$R_{C}$ for the both the outer and inner regions, and this result holds for
all five geometries. Although the equatorial geometries display enhanced
deficits of material in their inner regions, the net effect on the total
column density is much smaller (see the following subsection), only a factor
of $\sim 2$, so that the corresponding effect on the spectral energy
distributions is similarly modest (see Sec. 5).

Figure 5: The column density of the envelope along a ray defined by the polar
angle $\theta$. The envelope density is integrated from the planet radius
$R_{p}$ to the Hill radius $R_{H}$, using the fiducial values given in Table
1. Note that as $\theta\rightarrow\pi/2$, the isotropic, quasiequatorial, and
equatorial infall distributions have a large column density due to
concentration at the edge of the disk.  Figure 6: The column density of the
envelope. The top row corresponds to a viewing angle of $0$, from above. The
middle row is viewed from an angle of $\pi/4$, and the bottom row is viewed
from an angle of $\pi/2$, from the equatorial plane. The color at each point
represents the column density of the envelope along the ray through the
system. The contours are log-spaced evenly between $10^{-4}$ and $10^{-1}$
with a step of $1/2$.

###  3.4 Column Density

Due to the complexity of the orbit equation (Eq. 2), finding the column
density at an arbitrary polar parameter $\mu=\cos\theta$ is not analytically
tractable. However, we can numerically integrate the density (Eq. 7) along
rays with different lines of sight to the center, which we present in Fig. 5.
For the polar and quasipolar flow functions, the column density is relatively
flat as a function of viewing angle, since the behavior of the streamlines
counteracts the mass concentration. Even the isotropic case shows only modest
variations with $\theta$. In contrast, the equatorial and quasiequatorial flow
functions result in much more significant angular dependence, most notably a
deficit of column density along the rotational poles of the system.

It is interesting to note that in Fig. 5, the polar inflow column density is
essentially flat. This is a result of the fact that the infall tracks
concentrate material towards the equatorial plane, while that infall
distribution sets the matter to preferentially fall along the poles. These
factors nearly exactly cancel. Similarly, the column densities of the other
inflow functions become large at the disk plane, since material is
concentrated both from the equatorial inflow and the concentration of the
material.

In order to illustrate the density structure of the system, Fig. 6 shows
column density maps of the entire envelope along parallel lines of sight
through the system. The figure shows three viewing angles — the top row is
viewed from above with $\theta=0$, the middle row is viewed from midlatitudes
with $\theta=\pi/4$, and the bottom row is viewed from the equator with
$\theta=\pi/2$. This figure shows the three-dimensional structure of the
envelope for the different inflow distributions. In the remainder of this
section, we present a simplified analytical solution, which is found by
integrating through the equivalent spherical density distribution derived in
Sec. 3.3, which provides an approximate column density for the whole system.

Integrating through the density distribution, we evaluate

$\displaystyle N_{\rm col}=$
$\displaystyle\int\displaylimits_{0}^{\infty}\langle\rho\rangle{\rm
d}r=\int\displaylimits_{0}^{\infty}Cr^{-3/2}\mathcal{A}(u){\rm d}r$ (20)
$\displaystyle=$ $\displaystyle
CR_{C}^{-1/2}\int\displaylimits_{0}^{\infty}\frac{\sqrt{2u}}{u^{3/2}}{\rm
d}u\int\displaylimits_{\eta}^{1}\frac{f(x)}{[a^{2}+\mu_{0}^{2}]^{1/2}}{\rm
d}\mu_{0}\,.$

For the final equality, we have used the definition of the asphericity
function $\mathcal{A}(u)$. We can therefore switch the order of the
integration, integrating the cosine of the angle $\mu_{0}$ over the range
[0,1], while the integral in the radial variable $u$ develops a nonzero lower
limit. We define $u_{\rm in}(x)=1-\mu_{0}^{2}$ and write Eq. (20) as

$N_{\rm col}=C\sqrt{2R_{C}}\int\displaylimits_{0}^{1}f(\mu_{0}){\rm
d}\mu_{0}\int\displaylimits_{u_{\rm in}}^{\infty}\frac{{\rm
d}u}{u\sqrt{2u-1+\mu_{0}^{2}}}\,.$ (21)

The second integral can be evaluated and is equivalent to
$(\pi/2)\,(1-\mu_{0}^{2})^{-1/2}$. Therefore, the column density is

$\displaystyle N_{\rm col}=$
$\displaystyle\sqrt{2}\frac{\pi}{2}CR_{C}^{-1/2}\int\displaylimits_{0}^{1}\frac{f(\mu_{0}){\rm
d}\mu_{0}}{\sqrt{1-\mu_{0}^{2}}}$ (22) $\displaystyle\equiv$
$\displaystyle\sqrt{2}\frac{\pi}{2}CR_{C}^{-1/2}J\,.$

The second equality has defined the integral $J$. Now, we simply must evaluate
$J$ for the different inflow functions defined by Eq. (1). We find that

$J=\begin{dcases}3\pi/4&\text{polar},\\\ 2&\text{quasipolar},\\\
\pi/2&\text{isotropic},\\\ 4/\pi&\text{quasieq.},\\\
3\pi/8&\text{equatorial}.\end{dcases}$ (23)

Given this result, the estimated visual extinction for a standard opacity of
$\kappa_{V}=200$ cm2 g-1 is

$A_{V}=3.3J\Big{(}\frac{\dot{M}}{1\,M_{J}/{\rm
Myr}}\Big{)}\Big{(}\frac{M}{M_{J}}\Big{)}^{-2/3}\,.$ (24)

##  4 Circumplanetary Disk Properties

With the properties of the infalling envelope specified, we now derive the
properties of the circumplanetary disk that forms at the center of the
collapsing flow.

###  4.1 Magnetic Truncation

In general, the inner boundary of the disk will not extend all the way to the
planetary surface. Instead, the disk is truncated by the planetary magnetic
field. To leading order, the planetary magnetic field will have a dipole form

$\vec{B}(r,\mu)=B_{p}\Big{(}\frac{R_{p}}{r}\Big{)}^{3}[3\mu\hat{r}-\hat{z}]\,.$
(25)

The truncation radius $R_{X}$ is given by the radial point in the disk where
the ram pressure of the inward accretion flow is balanced by the outward
pressure from the magnetic field. This magnetic truncation radius has the form

$R_{X}=\omega\bigg{(}\frac{B_{p}^{4}R_{p}^{12}}{GM_{p}\dot{M}^{2}}\bigg{)}^{1/7}\,,$
(26)

where $\omega$ is a dimensionless constant of order unity (for derivations of
this equation, see Ghosh and Lamb 1978; Blandford and Payne 1982). Inserting
typical numbers, the truncation radius can be written as a scaled equation in
the form

$\displaystyle\frac{R_{X}}{R_{p}}\simeq 3.8$
$\displaystyle\Big{(}\frac{M}{M_{J}}\Big{)}^{-1/7}\Big{(}\frac{\dot{M}}{1\,M_{J}/{\rm
Myr}}\Big{)}^{-2/7}$ (27) $\displaystyle\times\Big{(}\frac{B_{p}}{500\text{
Gauss}}\Big{)}^{4/7}\Big{(}\frac{R_{p}}{10^{10}\text{ cm}}\Big{)}^{5/7}\,.$

###  4.2 Surface Density for Direct Infall

We first consider the surface density of the disk in the absence of viscous
mass transfer. In this limit, the surface density at a cylindrical radius $r$
gains mass at the rate at which incoming material hits the midplane of the
system. This rate is given by $2\rho v_{\theta}$, where the factor of $2$
arises because incoming material strikes both sides of the disk. Evaluating
this quantity at $\mu=0$ using Eqs. (2), (5), and (7), we find that

$\frac{{\rm d}\Sigma}{{\rm d}t}=2\rho v_{\theta}=\frac{\dot{M}}{4\pi
r^{2}}\frac{f_{i}(\mu_{0})}{\zeta\mu_{0}}=\frac{\dot{M}}{4\pi
r^{2}}\frac{f_{i}(\mu_{0})}{(\zeta^{2}-\zeta)^{1/2}}\,.$ (28)

The surface density built up over time is therefore given by the integral

$\Sigma_{0}(r)=\int\frac{\dot{M}}{4\pi
r^{2}}\frac{f_{i}(\mu_{0})}{(\zeta^{2}-\zeta)^{1/2}}{\rm d}t\,.$ (29)

Note that we can change variables by writing ${\rm d}M=\dot{M}{\rm d}t$. Since

$\zeta=\frac{R_{H}}{3r}=\frac{a}{3r}\Big{(}\frac{M}{3M_{\star}}\Big{)}^{1/3}\,,$
(30)

we can also write the surface density integral in the form

$\Sigma_{0}=\frac{243rM_{\star}}{4\pi
a^{3}}\int\displaylimits_{1}^{R_{C}/r}\frac{f_{i}(\mu_{0})\zeta^{2}}{(\zeta^{2}-\zeta)^{1/2}}{\rm
d}\zeta\,.$ (31)

Given that $r<R_{C}$ and $\mu=0$, we can use Eq. (2) to write
$\mu_{0}=\sqrt{1-1/\zeta}$.

For the isotropic case, $f_{i}(\mu_{0})=1$, so we can evaluate the integral in
Eq. (31) to find that

$\displaystyle\Sigma_{0}=\frac{3M}{16\pi R_{C}^{2}}\Big{[}$
$\displaystyle\frac{(2+3u)\sqrt{1-u}}{u}$ (32)
$\displaystyle+3u\tanh^{-1}\left(\sqrt{1-u}\right)\Big{]}\,.$

We have once again written the result in terms of $u=1/\zeta=r/R_{C}$.

For the quasipolar case, $f_{i}(\mu_{0})=2\mu_{0}=2\sqrt{(\zeta-1)/\zeta}$. As
a result, evaluating Eq. (31) yields

$\Sigma_{0}=\frac{3M}{4\pi R_{C}^{2}}\Big{[}\frac{1}{u}-u\Big{]}\,.$ (33)

For the polar case, $f_{i}(\mu_{0})=3\mu_{0}^{2}=3(\zeta-1)/\zeta$, so that
Eq. (31) evaluates to

$\displaystyle\Sigma_{0}=\frac{9M}{16\pi R_{C}^{2}}\Big{[}$
$\displaystyle\frac{(2-u)\sqrt{1-u}}{u}$ (34)
$\displaystyle-u\tanh^{-1}\left(\sqrt{1-u}\right)\Big{]}\,.$

For the quasiequatorial case,
$f_{i}(\mu_{0})=4\sqrt{1-\mu_{0}^{2}}/\pi=4\zeta^{-1/2}/\pi$, and we obtain

$\Sigma_{0}=\frac{2M}{\pi^{2}R_{C}^{2}}(1+2u)\sqrt{\frac{1-u}{u}}\,.$ (35)

Finally, the equatorial case gives
$f_{i}(\mu_{0})=3(1-\mu_{0}^{2})/2=3/2\zeta$, so we evaluate Eq. (31) to find

$\Sigma_{0}=\frac{9M}{8\pi
R_{C}^{2}}\left(\sqrt{1-u}+u\tanh^{-1}(\sqrt{1-u})\right)\,.$ (36)

As a consistency check, we can integrate the surface density over the entire
surface, that is, $\int\displaylimits_{0}^{1}2\pi R_{C}^{2}\Sigma_{0}u{\rm
d}u$. For all of these functions, we obtain the total mass $M$ which has come
through the Hill sphere, as expected. Summarizing, we write the surface
density distribution as

$\Sigma_{0}=\frac{3M}{4\pi R_{C}^{2}}\begin{dcases}\begin{aligned}
\frac{3}{4}\Big{[}&\frac{(2-u)\sqrt{1-u}}{u}\\\
-&u\tanh^{-1}\left(\sqrt{1-u}\right)\Big{]}\end{aligned}&\text{polar,}\\\
\frac{1-u^{2}}{u}&\text{quasipolar,}\\\ \begin{aligned}
\frac{1}{4u}\big{[}&(2+3u)\sqrt{1-u}\\\
+&3u^{2}\tanh^{-1}(\sqrt{1-u})\big{]}\end{aligned}&\text{isotropic,}\\\
\frac{8}{3\pi}\frac{(1+2u)(1-u)^{1/2}}{\sqrt{u}}&\text{quasieq.,}\\\
\begin{aligned} \frac{3}{2}\big{[}&\sqrt{1-u}\\\
+&u\tanh^{-1}(\sqrt{1-u})\big{]}\end{aligned}&\text{equatorial.}\end{dcases}$
(37)

In the limit where $u\rightarrow 0$, the isotropic and both polar solutions
approach $\Sigma_{0}\propto u^{-1}$, the quasiequatorial solution approaches
$\Sigma_{0}\propto u^{-1/2}$, and the equatorial solution approaches a
constant, which grows even more slowly. Fig. 7 shows the behavior of these
surface density functions over the radial extent of the disk.

Figure 7: Surface density distributions for forming circumplanetary disks for
the case of direct infall only. The curves show the results for varying input
flow functions, from polar (top, red) to equatorial (bottom, purple). All of
the surface density distributions have $\Sigma_{0}\rightarrow 0$ as
$r\rightarrow R_{C}$, because no material impacts the equatorial plane beyond
the outer edge of the disk. The inner edge is at $r=R_{X}$, the truncation
radius.

###  4.3 Steady-State Surface Density Distributions

Now, we include viscous accretion as a modification to the surface density.
The equation of motion for a viscously-evolving disk has the form

$\displaystyle\frac{\partial\Sigma}{\partial t}=$
$\displaystyle\frac{3}{r}\frac{\partial}{\partial
r}\bigg{[}r^{1/2}\frac{\partial}{\partial
r}\Big{(}r^{1/2}\nu\Sigma\Big{)}\bigg{]}$ (38)
$\displaystyle+\frac{\dot{M}}{4\pi
r^{2}}\frac{f_{i}(\mu_{0})}{(\zeta^{2}-\zeta)^{1/2}}\,.$

Eq. (38) includes the infall source term, and $\nu=\alpha c_{s}^{2}/\Omega$ is
the kinematic viscosity. For a steady-state solution, the time derivatives
vanish. Defining $S\equiv 6\pi\nu\Sigma/\dot{M}$, Eq. (38) reduces to the
steady-state form

$2\frac{\partial}{\partial u}\bigg{[}u^{1/2}\frac{\partial}{\partial
u}\Big{(}u^{1/2}S\Big{)}\bigg{]}+\frac{f_{i}(\mu_{0})}{(1-u)^{1/2}}=0\,,$ (39)

where $u=r/R_{C}=1/\zeta$ as usual. Note that the mass accretion rate through
the disk is related to the function $S$ by

$\dot{M}_{\rm disk}=\dot{M}u^{1/2}\frac{\partial}{\partial
u}\left(u^{1/2}S\right)\,,$ (40)

where $\dot{M}$ is the net rate at which the planet/disk systems gains mass.
For each inflow function, we can integrate Eq. (39) to find the solution for
the reduced surface density $S(u)$, with integration constants $K_{1}$ and
$K_{2}$, that are chosen to satisfy the boundary conditions.

For isotropic inflow, we find that

$S=\frac{1}{\sqrt{u}}\left[\sin^{-1}\sqrt{u}+\sqrt{u(1-u)}\right]+2K_{1}+\frac{K_{2}}{\sqrt{u}}\,.$
(41)

For the quasipolar inflow function, we find that

$S=2\left[1-\frac{u}{3}\right]+2K_{1}+\frac{K_{2}}{\sqrt{u}}\,.$ (42)

For the polar case, the reduced density profile has the form

$\displaystyle S=\frac{1}{4}\Big{[}$
$\displaystyle(5-2u)\sqrt{1-u}+\frac{3}{\sqrt{u}}\sin^{-1}\sqrt{u}\Big{]}$
(43) $\displaystyle+2K_{1}+\frac{K_{2}}{\sqrt{u}}\,.$

For the quasiequatorial case, the profile has the form

$\displaystyle S=\frac{4}{\pi}\frac{1}{\sqrt{u}}\Big{[}$
$\displaystyle\sqrt{u}\cos^{-1}\sqrt{u}+\frac{4}{3}-\frac{1}{3}(4-u)\sqrt{1-u}\Big{]}$
(44) $\displaystyle+2K_{1}+\frac{K_{2}}{\sqrt{u}}\,.$

Finally, the equatorial surface density profile takes the form of

$S=\frac{1}{8}\left[(7+2u)\sqrt{1-u}+\frac{9}{\sqrt{u}}\sin^{-1}\sqrt{u}\right]+2K_{1}+\frac{K_{2}}{\sqrt{u}}\,.$
(45)

Ignoring the $K_{1}$ and $K_{2}$ terms, the steady-state surface density can
be written in the form

$\Sigma=\frac{\dot{M}}{6\pi\nu}\begin{dcases}\begin{aligned}
\frac{1}{4\sqrt{u}}\big{[}&(5-2u)\sqrt{u(1-u)}\\\
&+3\sin^{-1}\sqrt{u}\big{]}\end{aligned}&\text{polar,}\\\
2(1-u/3)&\text{quasipolar,}\\\
\frac{1}{\sqrt{u}}\big{[}\sqrt{u(1-u)}+\sin^{-1}\sqrt{u}\big{]}&\text{isotropic,}\\\
\begin{aligned} \frac{4}{3\pi\sqrt{u}}\big{[}&3\sqrt{u}\cos^{-1}\sqrt{u}\\\
&+4-(4-u)\sqrt{1-u}\big{]}\end{aligned}&\text{quasieq.,}\\\ \begin{aligned}
\frac{1}{8\sqrt{u}}\big{[}&(7+2u)\sqrt{u(1-u)}\\\
&+9\sin^{-1}\sqrt{u}\big{]}\end{aligned}&\text{equatorial.}\end{dcases}$ (46)

These profiles are shown in Fig. 8, where we assume an $\alpha$-viscosity
profile (Shakura and Sunyaev, 1973).

We can then evaluate Eq. (40) to find the limits of the accretion rate. As
these limits are equivalent for all of the inflow functions, we will discuss
isotropic inflow and then generalize to the other functional forms. For
isotropic inflow,

$\dot{M}_{\rm disk}=\dot{M}\big{[}\sqrt{1-u}+K_{1}\big{]}\,.$ (47)

In the limit $K_{1}\rightarrow 0$, the disk accretion rate is thus equal to
the mass inflow into the entire system (as $u\rightarrow 0$). Since some
fraction of the inflowing mass is accreted directly onto the planet, this
provides a boundary condition for $K_{1}$, to ensure that the correct fraction
is accreted through the disk. As $u\rightarrow 1$, the mass accretion rate
vanishes, since there is no material beyond $u=1$ to be brought inward. As
mentioned, the other profiles behave identically at these limits with only
minor differences in the leading terms.

It is interesting to note that the surface density due to direct accretion
(Fig. 7) shows far more significant qualitative differences among the
different flow functions when compared to the surface density for a steady-
state viscous disk (Fig. 8). As long as the viscosity is large enough, so that
the system can reach steady-state, disk evolution smooths out the surface
density into a nearly universal form. The viscous time is given by
$t_{vis}\sim R_{C}^{2}/\nu$, which must be shorter than the time over which
the disk mass changes from infall $t_{inf}\sim M/{\dot{M}}$. If we evaluate
the viscous accretion time for a typical system (a planet at $a$ = 5 au in a
minimum-mass solar nebula), then $t_{vis}\approx$ (16 yr)/$\alpha$, which will
be much shorter than $t_{inf}\sim 1$ Myr, as long as $\alpha>10^{-4}$ (see
Adams and Batygin 2022 for a more detailed discussion). We thus expect the
circumplanetary disks to reach steady-state conditions.

Figure 8: Steady-state surface density in a viscous disk, for the different
infall functions. While the surface density stops at the edge of the disk,
$r=R_{C}$, practically the density should continue outside of this limit. This
is due to the fact that angular momentum conservation will necessarily expand
the disk as viscosity operates, moving material beyond $R_{C}$. The inner edge
is at $r=R_{X}$, the truncation radius.

###  4.4 Direct Accretion onto the Planet

In the core accretion paradigm, a certain fraction of the infalling matter is
accreted directly onto the planetary surface. In this section, we calculate
the fraction of material that falls onto the surface, and the planetary mass
that would result in the absence of disk accretion. This mass scale is
generally much smaller than the expected planetary masses, and this
discrepancy underscores the importance of disk accretion (most likely due to
viscosity).

The rate of direct infall through a sphere of radius $r<R_{H}$ is given by

$\dot{M}_{\rm direct}=4\pi r^{2}\int\displaylimits_{0}^{1}\rho|v_{r}|{\rm
d}\mu=\dot{M}\int\displaylimits_{\eta}^{1}f_{i}(\mu_{0}){\rm d}\mu_{0}\,.$
(48)

This result is derived using the definition of the density (Eq. 7), along with
the integral limit $\eta$ given by Eq. (9).

We first consider the rate of direct infall onto the planet’s surface
directly, defining a fraction $f_{p}$ such that $\dot{M}_{\rm
direct}\equiv\dot{M}f_{p}$, with $f_{p}$ given by the integral in Eq. (48). By
noting that $\eta=\sqrt{1-u}$ (from Eq. 9), we can determine that
$\eta=\sqrt{1-R_{p}/R_{C}}$. Evaluating the $f_{p}$ integral with the inflow
functions given by Eq. (1), we find that

$f_{p}=\begin{dcases}1-(1-u_{p})^{3/2}&\text{polar,}\\\
u_{p}&\text{quasipolar,}\\\ 1-\sqrt{1-u_{p}}&\text{isotropic,}\\\
\begin{aligned} &1-(2/\pi)\Big{[}\sqrt{u_{p}(1-u_{p})}\\\
&+\arcsin{\sqrt{1-u_{p}}}\Big{]}\end{aligned}&\text{quasieq.,}\\\
1-\sqrt{1-u_{p}}(1+u_{p}/2)&\text{equatorial.}\end{dcases}$ (49)

Note that $u_{p}=R_{p}/R_{C}$. In addition, $f_{d}$, the fraction of material
that flows directly onto the disk, is given by $f_{d}=1-f_{p}$. Working in the
limit $R_{C}\gg R_{p}$, we can determine the leading order expressions for the
direct infall fractions for the flow functions of interest, i.e.,

$f_{p}\simeq\begin{dcases}3u_{p}/2&\text{polar,}\\\
u_{p}&\text{quasipolar,}\\\ u_{p}/2&\text{isotropic,}\\\
(4/3\pi)u_{p}^{3/2}&\text{quasieq.,}\\\
(3/8)u_{p}^{2}&\text{equatorial.}\end{dcases}$ (50)

Assuming that the inflow function is consistent over the planet’s lifetime,
the total mass of the planet is given by

$M_{p}=\int\displaylimits_{0}^{t}f_{p}\dot{M}\,{\rm
d}t=M_{0}+\int\displaylimits_{M_{0}}^{M}f_{p}\,{\rm d}M\,.$ (51)

In Eq. (51), $M_{0}$ is an initial mass value, defined to be the mass where
the centrifugal radius is equal to the planetary radius, i.e.,

$M_{0}=81M_{\star}\left(\frac{R_{p}}{a}\right)^{3}\,.$ (52)

However, the runaway gas accretion phase of giant planet formation is
generally understood to begin when the planetary core mass is $M_{p}\simeq 20$
${\rm M}_{\oplus}$. Since the disk-forming mass given by Eq. (52) may be
smaller than $20$ ${\rm M}_{\oplus}$, caution must be taken in determining the
initial mass. However, in either case we will assume that the final mass is
large compared to the initial mass, and so this correction is small.

Regardless, we will use the fact that $R_{C}\gg R_{p}$ to replace the $f_{p}$
in Eq. (51) with the approximate expressions in Eq. (50). We can then use Eq.
(52) and the fact that $R_{C}=(a^{3}M/81M_{\star})^{1/3}$ to show that

$\frac{R_{p}}{R_{C}}=\left(\frac{M_{0}}{M}\right)^{1/3}\,.$ (53)

We use this relation to replace the $R_{p}/R_{C}$ in Eq. (50). Integrating and
assuming that $M_{0}\ll M_{p}$, we find that

$M_{p}\simeq\begin{dcases}\frac{9}{4}(M_{0}M^{2})^{1/3}&\text{polar,}\\\
\frac{3}{2}(M_{0}M^{2})^{1/3}&\text{quasipolar,}\\\
\frac{3}{4}(M_{0}M^{2})^{1/3}&\text{isotropic,}\\\
\frac{8}{3\pi}(M_{0}M)^{1/2}&\text{quasieq.,}\\\
\frac{9}{8}(M_{0}^{2}M)^{1/3}&\text{equatorial.}\end{dcases}$ (54)

Keep in mind that $M$ is the total mass of the planet–disk system, whereas
$M_{p}$ is the mass of the planet. It is important to note, however, that this
accounts only for accretion onto the planet’s surface directly, and does not
consider the additional effect of magnetic fields, which will expand the
region where infalling matter is funneled directly onto the planet (Adams and
Batygin, 2022). Incorporation of the magnetic fields, which becomes
complicated, is left for future work.

##  5 Radiation Signatures

In this section, we find the luminosity of the various components of the
planet–disk–envelope system. As outlined below, determining the total system
luminosity is complicated by the somewhat-large number of different
contributions. For the expected case where the disk viscosity is sufficiently
large, however, most of the incoming material falls onto the disk surface but
eventually accretes onto the central planet. The luminosity is thus dominated
by the viscous contribution, along with any intrinsic luminosity within the
planet itself. With the luminosity specified, we can estimate the radiative
appearance of these forming planets. Here we present an approximate treatment,
but leave the full two-dimensional radiative transfer problem for a
forthcoming paper.

###  5.1 Planet Luminosity

First, let us consider the luminosity of the planet. We will assume that the
planet is in thermal equilibrium, with a well-defined outer radius. As a
result, the planet reemits all the power due to accretion that is delivered to
its surface (where we neglect the planetary winds, which could carry some of
the power). The benchmark power delivered to the entire system is the energy
generated by all the material falling onto the planet’s surface and
dissipating its kinetic energy as if it fell from infinity. The energy scale
is determined by the gravitational potential of the planet, so that the
benchmark power is given by

$L_{0}\equiv\frac{GM\dot{M}}{R_{p}}\,.$ (55)

A certain fraction of the incoming material $f_{p}$ falls onto the planet,
while a fraction $f_{d}$ falls onto the disk. The planet–falling fraction
$f_{p}$ is given by Eq. (49), while the disk-falling fraction is
$f_{d}=1-f_{p}$. Note that aside from the values of the fractions
$(f_{p},f_{d})$, the following results are general for all five inflow
functions.

The planet is assumed to be corotating with the disk at the magnetic
truncation radius $R_{X}$, where the Keplerian rotation rate is given by

$\Omega_{p}^{2}=\frac{GM_{p}}{R_{X}^{3}}\,.$ (56)

As a result, the matter on the planetary surface has an azimuthal rotation
velocity given by

$v_{\phi}=R_{p}\Omega_{p}\sin\theta\,.$ (57)

In Eq. (57), $\theta$ is the colatitude on the planetary surface. For a given
mass-inflow rate, the rate of change of kinetic energy stored in rotation has
the form

$\frac{{\rm d}E_{\rm rot}}{{\rm
d}t}=\frac{1}{2}\frac{R_{p}^{3}}{R_{X}^{3}}\langle\sin^{2}\theta\rangle
L_{0}\,.$ (58)

Assuming that the infalling material is well-mixed on the planet’s surface,
the average of the geometric factor $\langle\sin^{2}\theta\rangle=2/3$, so
that the luminosity produced by the planet as a result of direct infall is

$L_{p}^{\rm dir}=f_{p}\Big{(}1-\frac{R_{p}^{3}}{3R_{X}^{3}}\Big{)}L_{0}\,.$
(59)

The luminosity produced via accretion from the disk onto the planet is
similar, with a small correction. Specifically, matter no longer falls from
infinity, but from the magnetic truncation radius $R_{X}$. As a result, the
planet luminosity from material accreted from the disk is given by

$L_{p}^{\rm
acc}=f_{d}\Big{(}1-\frac{R_{p}^{3}}{3R_{X}^{3}}\Big{)}\Big{(}1-\frac{R_{p}}{R_{X}}\Big{)}L_{0}\,.$
(60)

Adding Eqs. (59) and (60), we find that the total luminosity of the planet is

$L_{p}=L_{0}\Big{(}1-\frac{R_{p}^{3}}{3R_{X}^{3}}\Big{)}\Big{(}1-f_{d}\frac{R_{p}}{R_{X}}\Big{)}+L_{\rm
int}\,,$ (61)

where $L_{\rm int}$ is the intrinsic luminosity of the planet, including both
power from gravitational contraction and radiogenic heating in the core. The
resulting planetary temperature is given by

$\sigma T_{p}^{4}=\frac{L_{0}}{4\pi
R_{p}^{2}}\Big{(}1-\frac{R_{p}^{3}}{3R_{X}^{3}}\Big{)}\Big{(}1-f_{d}\frac{R_{p}}{R_{X}}\Big{)}+\frac{L_{\rm
int}}{4\pi R_{p}^{2}}\,.$ (62)

In general (and specifically under the approximations of this paper), the
planet’s internal luminosity is small compared to the other two terms, and so
that $L_{\rm int}$ can be dismissed from Eqs. (61) and (62).

Assuming that the planet is a blackbody, the wavelength-dependent planet
luminosity is given by

$L_{p,\nu}=4\pi R_{p}^{2}\,\pi B_{\nu}(T_{p})\,,$ (63)

and the corresponding monochromatic flux takes the form

$F_{p,\nu}(r,\mu)=\frac{L_{p,\nu}}{4\pi
r^{2}}G_{pd}(\mu)\exp[-\tau_{\nu}(r,\mu)]\,.$ (64)

In Eq. (64), $\tau_{\nu}$ is the optical depth of the envelope between the
point $(r,\mu)$ and the planet, which depends on the opacity and the column
density. The function $G_{pd}(\mu)$ accounts for the possible shadowing of the
disk by the planet, but is set to unity in the present treatment. Note that
while the column density of the envelope is exactly specified by the solutions
of the previous section, the density structure is complicated and the column
density must be evaluated numerically.

###  5.2 Disk Luminosity

Next, we calculate the disk luminosity. The viscosity is expected to be large
enough that most of the material landing on the disk surface will eventually
accrete onto the planet (see also Zhu 2015; Zhu et al. 2018). As a result, the
total disk luminosity is $L_{d}=f_{d}L_{0}$, but some fraction of this power
is already accounted for by the material that falls from the magnetic
truncation radius to the planetary surface (specifically $L_{p}^{\rm acc}$).
In addition, half of the energy remaining in the disk is stored in rotation
due to the virial theorem, so the total disk luminosity is given by

$L_{d}=f_{d}L_{0}\frac{R_{p}}{2R_{X}}\,.$ (65)

With the disk luminosity specified, we need to estimate the surface
temperature of the disk. For both a steady-state accretion disk and a flat,
passively heated disk, the temperature distribution is a power-law of the form

$T_{d}(r)=T_{X}\Big{(}\frac{R_{X}}{r}\Big{)}^{3/4}\,.$ (66)

Since the total luminosity of the disk is $L_{d}$, conservation of energy
implies that

$L_{d}=2\pi\int\displaylimits_{R_{X}}^{R_{C}}\sigma T_{d}^{4}(r)r\,{\rm
d}r\,.$ (67)

The integral in Eq. (67) can be solved with Eqs. (65) and (66) to find that

$\sigma T_{X}^{4}=f_{d}\frac{GM\dot{M}}{8\pi
R_{X}^{3}}\Big{(}1-\frac{R_{X}}{R_{C}}\Big{)}^{-1}\,.$ (68)

Note that we have implicitly assumed that the majority of the disk luminosity
is produced by viscosity, rather than by the flow of material directly onto
the disk. This approximation holds because the energy released by viscosity is
the gravitational potential energy of falling from the disk onto the planet.
For comparison, the energy released by the material striking the disk is the
gravitational potential energy that falls onto the disk from infinity. The
disk radius is much larger than the planet radius, so that the energy released
by viscosity dominates over the energy released by the direct mass inflow.
Moreover, the luminosity due to shocks on the disk surface and mixing of the
newly added material has been calculated previously (Adams and Batygin 2022;
see also Cassen and Moosman 1981; Adams and Shu 1986).

The total wavelength-dependent luminosity of the disk is given by

$L_{d,\nu}=2(2\cos\theta)\int\displaylimits_{R_{X}}^{R_{C}}\pi
B_{\nu}[T_{d}(r)]2\pi r\,{\rm d}r\,.$ (69)

We assume that the disk is optically thick to its internal radiation. The
factor of $2$ arises from the two sides of the disk, and the factor of
$(2\cos\theta)$ arises from a reduction of the effective area due to the
viewing angle $\theta$, which must be normalized. The corresponding
monochromatic flux takes the form

$F_{d,\nu}(r,\mu)=\frac{1}{2}\frac{L_{d,\nu}}{4\pi
r^{2}}G_{dp}(\mu)\exp(\tau_{\nu}(r,\mu))\,,$ (70)

where $G_{dp}(\mu)$ accounts for the shadowing of the disk by the planet, and
$\tau_{\nu}$ is the optical depth of the envelope between the observation
point and the disk. Note that there is a factor of $1/2$, which accounts for
the fact that we assume that the disk is optically thick to itself, and so
always shadows one of its faces. In addition, the use of a simple $\tau_{\nu}$
term for the optical depth represents an approximation — because the disk is
physically extended, the optical depth should be computed between every point
on the disk surface and the observer, with the radiation from each surface
element reduced by an appropriate factor. This integral is well-defined by our
density solutions but is complicated to solve analytically, and is left for
numerical evaluation in our forthcoming paper.

###  5.3 Envelope Emission

With the luminosity contributions specified, we estimate the spectral energy
distribution for the forming planet, including the infalling envelope, the
disk, and the planet itself. The envelope has no intrinsic luminosity, but it
intercepts radiation from the central planet and disk and reradiates the
energy at longer wavelengths. The envelope is (mostly) optically thick to the
infrared radiation emitted by the planet and disk. In contrast, the envelope
temperatures are lower and the reemitted radiation has longer wavelengths, so
the envelope is generally optically thin to its own radiation. Since the
envelope is rotating and not spherically symmetric, the resulting radiative
transfer problem is two-dimensional. Leaving a full exploration of this
problem for future work, we here present an approximate treatment.

Toward this end, we use the equivalent spherical envelope approximation from
Sec. 3.3. Since the envelope is optically thin to its own radiation, the most
important property is the total column density (or optical depth) and not the
specific geometry. Due to conservation of energy, the total effective envelope
luminosity $L_{e}$ is determined by the fraction of the central source
radiation that is attenuated by the infalling envelope, i.e.,

$L_{e}=L_{p}+L_{d}-\int\displaylimits_{0}^{\infty}(L_{p,\nu}+L_{d,\nu})\exp(-\kappa_{\nu}N_{\rm
col})\,{\rm d}\nu\,.$ (71)

Note that in Eq. (71), we have not only assumed spherical symmetry of the
envelope, but placed the disk at the center of the system. In reality, of
course, the disk is physically extended across the equatorial plane to $1/3$
of the envelope radius, but most of the radiation is emitted from the inner
regions of the disk. In addition, the radiation emitted from the outer disk
regions has a longer wavelength and is less susceptible to being absorbed.

The effective luminosity $L_{e}$ of the envelope is reemitted at longer
wavelengths. The total luminosity is given by the integral over the envelope
volume and takes the form

$L_{e}=16\pi\int\displaylimits_{R_{p}}^{R_{H}}\rho\kappa_{P}\sigma
T_{e}^{4}r^{2}\,{\rm d}r\,,$ (72)

where $\kappa_{P}$ is the Planck mean opacity. One factor of $4\pi$ comes from
the surface area term, and another comes from integration over solid angle
(partially cancelled by the definition of the Stefan-Boltzmann constant
$\sigma$). Because the temperature range of the envelope is limited (ranging
from $T_{p}\sim 1000$ K at the planetary surface to $T_{d}\sim 100$ K at the
Hill radius), we can assume a simple power law for the opacity,

$\kappa_{\nu}=\kappa_{0}\Big{(}\frac{\nu}{\nu_{0}}\Big{)}^{\eta}\,.$ (73)

The index is expected to fall in the range $1\leq\eta\leq 2$ (see Draine and
Lee 1984 and subsequent references). With this form of the frequency dependent
opacity, the Planck mean opacity can be analytically evaluated and takes the
form

$\kappa_{P}\equiv
b_{\kappa}T^{\eta}=\kappa_{0}\Big{(}\frac{k_{B}}{\nu_{0}h}\Big{)}^{\eta}\frac{\Gamma(4+\eta)\zeta_{R}(4+\eta)}{6\zeta_{R}(4)}\,T^{\eta}\,,$
(74)

where $\zeta_{R}$ is the Riemann zeta function and $\Gamma$ is the gamma
function (Abramowitz and Stegun, 1972). Solving the radiative transfer
equations with spherical symmetry (Adams and Shu, 1985), and adopting an index
$\eta=1$ (see also Semenov et al. 2003), the temperature distribution in the
envelope takes the form

$T_{e}(r)=T_{C}\left(\frac{r}{R_{C}}\right)^{-2/5}\,.$ (75)

With this temperature profile and the power-law form for the opacity, the
luminosity integral (Eq. 72) can be evaluated to yield

$L_{e}=16\pi R_{C}^{2}\sigma b_{\kappa}T_{C}^{5}N_{\rm col}\,,$ (76)

where $N_{\rm col}$ is the total column density in the spherically-equivalent
limit, as defined by Eq. (20). However, that these are not exactly equivalent,
since Eq. (20) defines the column density to be the integral from $0$ to
$\infty$, while Eq. (72) integrates from the planet radius $R_{p}$ to the Hill
radius $R_{H}$. We use the latter result here.444Note that we are calculating
the SED at the Hill radius, although additional material in the background
circumstellar disk, beyond $R_{H}$, can process the radiation further. This
correction is situationally dependent, but can be included in future
applications.

Using the amended specification of the column density $N_{\rm col}$, we can
solve Eq. (76) for the fiducial temperature $T_{C}$ to obtain

$T_{C}=\bigg{(}\frac{L_{e}}{16\pi R_{C}^{2}\sigma b_{\kappa}N_{\rm
col}}\bigg{)}^{1/5}\,.$ (77)

With the fiducial temperature determined, we now have the complete (spherical)
temperature profile. The total monochromatic luminosity for the envelope is
then given by

$L_{e,\nu}=4\pi\int\displaylimits_{R_{X}}^{R_{H}}\langle\rho\rangle(r)\kappa_{\nu}4\pi
B_{\nu}[T_{e}(r)]r^{2}\,{\rm d}r\,.$ (78)

where $\langle\rho\rangle$ is the density field of the equivalent spherical
infall region, given by Sec. 3.3. The monochromatic flux is given by

$F_{e,\nu}=\frac{L_{e,\nu}}{4\pi r^{2}}\,.$ (79)

Note that there is no attenuation term, since the envelope is assumed to be
optically thin to its own radiation.

###  5.4 Spectral Energy Distributions

Using the above results, we can now estimate the spectral energy distributions
(SEDs) for forming planets. Keep in mind that we have assumed that (i) the
planet and the disk are both at the center of the envelope and (ii) the
envelope is spherically symmetric, using the concept of an equivalent
spherical envelope.

The SED is defined to be $\nu L_{\nu}$, where $\nu$ is the frequency. The SEDs
for the representative values shown in Table 1 are shown in Fig. 9. It is
worth noting that the disk luminosity (Eq. 69) depends on the viewing angle
$\theta$. For the purposes of this figure, we assume that $\theta=0$, so that
the system is seen from the rotational pole. Variation in the viewing angle is
generally an $\mathcal{O}(1)$ effect on the disk SED, only becoming
significant in the limit $\theta\rightarrow\pi/2$. The viewing angle has no
effect on the planet or envelope contribution to the SEDs in the spherically-
symmetric approximation used here.

The spectral energy distributions presented in Fig. 9 show several interesting
trends. The most significant differences appear in the envelope, where the
distinctions in the envelope density distribution (shown in Fig. 3) dominate.
For the parameter values used here, the disk emission contribution is
relatively small compared to that of the planet and the envelope. This finding
is in apparent contrast to Szulágyi et al. (2019), who found that the disk is
brighter than the planet at near infrared wavelengths. This difference is
likely the result different choices of the inner boundary condition,
especially the magnetic truncation radius, which determines the amount of
energy dissipated as incoming material falls from the inner disk edge onto the
planetary surface.

The luminosity of the planet and the disk are primarily set by the delivered
accretion power, which depends on the mass inflow rate $\dot{M}$ (which is the
same for all cases shown) and the fraction of matter delivered to the
planetary surface $f_{p}$ or equivalently, that delivered to the disk
$f_{d}=1-f_{p}$ (see Eq. 61). For the different inflow geometries, differences
in the fraction $f_{d}$ are order $\sim u_{p}\sim 10^{-2}$, so that this
luminosity does not vary greatly. On the other hand, the differences in column
density are of order unity (see Eq. 23). For these parameter values, the total
optical depth of the envelope and the corresponding envelope luminosity are
modestly large, so that the envelope luminosity is comparable to the
unattenuated luminosity of the planet. As a result, the relative contributions
of the planet and envelope vary with the infall geometry, as shown in the Fig.
9. The envelope contribution is more significant for the polar cases,
primarily due to increased optical depth.

Figure 9: The spectral energy distributions of the planet–disk–envelope
system, operating under the assumptions that (i) the envelope is spherically
symmetric, (ii) the disk is a point source at the center of the envelope, and
(iii) the envelope is optically thin to its own radiation. SEDs are shown for
all 5 inflow functions. The planet and system parameter values are given in
Table 1. The planet, disk, envelope, and total SEDs are all shown. The viewing
angle is assumed to be normal to the disk surface.

##  6 Conclusion

In this paper, we have developed an analytical model for the late stages of
planet formation according to the core accretion paradigm, expanding on the
work of Adams and Batygin (2022). Those authors found theoretical results for
the structure and radiative signatures of the circumplanetary disk and
infalling envelope associated with a forming giant planet, assuming that the
incoming mass flow is distributed isotropically. This paper extends those
results to include different inflow distributions at the Hill radius, which
marks the outer boundary of the system. These different boundary conditions
are indicated by different numerical simulations of forming planets (compare
Ayliffe and Bate 2009, Lambrechts and Lega 2017, and Szulágyi and Mordasini
2017, which show qualitative and quantitative differences in the inflow
distributions). By constructing theoretical models of this later stage of
planet formation, we obtain signatures of the distribution of infalling matter
that may be detectable in future observations of giant protoplanets.
Importantly, our analytical construction means that the solutions are
applicable over a wide range of parameter space, independent of the specific
canonical values used herein.

The results of this paper show that the varying infall geometries (envelope
geometries) lead to significant qualitative differences in the density
structure of the infalling envelope, where the differences are greatest within
the centrifugal radius (Fig. 5). The equidensity contours become prolate for
infall along the polar directions and oblate for infall along the equatorial
directions. The corresponding total column densities through the envelope show
variations of order unity (Sec. 4.4). The disk surface density distribution
resulting from direct accretion shows corresponding distinctions (Fig. 7),
with steeper profiles for infall along the poles. In the presence of
sufficiently large viscosity, however, these differences are smoothed out in
the steady-state solutions of the surface density functions (Fig. 8). As a
result, the steady-state disk surface density approaches a nearly universal
form as long as viscous accretion can keep up with the incoming mass flow.

This paper also presents a preliminary assessment of the radiative appearance
of forming planets for the different infall geometries under consideration. We
calculate the spectral energy distributions using the ansatz of an equivalent
spherical envelope, and find relatively modest observational differences (see
Fig. 9). In general, the column density of the envelope is larger for flow
geometries concentrated along the poles of the system (rather than the
equatorial regions). As a result, the polar cases show greater luminosity from
the infalling envelope, and less (unattenuated) radiation from the central
planet and disk. As a result, the SEDs for the equatorial cases have a double-
peaked morphology, while the polar cases vary more smoothly with frequency.

For the canonical set of planetary parameters depicted in Fig. 9, the total
optical depth of the envelope is of order unity, so that differences in the
density structure have a modest effect on the SEDs. For systems with smaller
optical depth (e.g., from smaller infall rate ${\dot{M}}$), comparatively
little radiation is absorbed by the envelope, so the resulting SEDs will be
similar for all geometries. For systems with larger optical depth, however,
the dependence on infall geometry becomes more important. More specifically,
with greater optical depth, the envelope will absorb and reradiate a greater
fraction of the central source luminosity and will display a larger effective
envelope luminosity. This reradiated emission depends on the density
structure, and the geometry will be especially important for studying resolved
systems. Notice also that for sufficiently large optical depths, the envelope
will no longer be optically thin to its internal radiation. As a result, a
full treatment of the envelope emission requires detailed radiative transfer
calculations in two spatial dimensions, and this task is left for future work.

The treatment of this paper makes a number of simplifying assumptions in order
to obtain analytic results. For example, the trajectories of the incoming
material are calculated in the pressure-free limit. The corrections due to
pressure forces are expected to be relatively modest under conditions of
efficient cooling (again see the analysis of Adams and Batygin 2022), but can
be substantial if the effective equation of state is sufficiently stiff. In
addition, we take ${\dot{M}}$ to be the net infall rate, thereby ignoring
material that enters the Hill sphere and immediately flows back out. This
extraneous material will not affect the forming planet itself, or the disk,
but its presence adds an additional component to the envelope, which in turn
affects the spectral appearance of the object. Finally, we have taken the
incoming material to have azimuthal symmetry. The outer boundary conditions
can break this symmetry through a variety of effects acting on the background
circumstellar disk, including gap opening by the planet, headwinds due to non-
Keplerian rotation, and spiral density waves. These complications will
generally affect the structure of the envelope more than the disk or planet,
but nonetheless should be addressed in future work.

A number of additional avenues for future research should be explored. This
present work does not include line emission, which may be especially important
when arising from the accretion shock on the planetary surface. This
radiation, especially the H$\alpha$ line, provides a metric for the system
luminosity, as well as a potentially powerful diagnostic of the shock
conditions on the planetary surface (e.g., Marleau et al. 2022, 2023). Another
issue that must be explored further is the effects of magnetic fields. Young
and forming planets are expected to have surface fields in the range 100 –
1000 G (e.g., Yadav and Thorngren 2017), large enough to affect the geometry
of the accretion flow onto the planet. In particular, the fields can enhance
the cross-section for directly capturing incoming material and thereby alter
the luminosity estimates. The structure of the circumplanetary disk does not
only affect mass flow onto the planet and the spectral appearance of the
system, but can influence the formation of satellites (Batygin and Morbidelli,
2020). Finally, when the growing planet becomes large enough to clear a gap in
the circumstellar disk flow into the Hill sphere can take place along
particular directions (see, e.g., Szulágyi et al. 2022), which will provide
even more complicated boundary conditions for the incoming material.

In this paper, we have considered the forming planet with mass $M_{p}$ (or
system mass $M$) to have a given mass infall rate ${\dot{M}}$, magnetic field
strength $B_{p}$, and semi-major axis $a$. For these values, we calculated the
solutions for the density and velocity fields of the infall, the column
density, the disk surface density, and so on, for each inflow geometry. In
general, however, the infall rate, the magnetic field strength, and even the
semi-major axis can change with time or — equivalently — as the planet grows
in mass. Since the crossing time of the Hill sphere is on the order of $\sim
1$ yr, much less than the expected evolution time of order $\sim 1$ Myr, our
approach using instantaneous values of the system variables is valid. We can
then incorporate evolutionary scenarios into this framework by specifying how
the system properties vary with mass, although the relevant functions
${\dot{M}}(M_{p})$, $B_{p}(M_{p})$, and $a(M_{p})$ are not yet fully
understood. In addition, the flow geometry can change as the planet grows. For
example, for a relatively small planet the Hill sphere will be safely smaller
than the scale height of the circumstellar disk, so that an isotropic or polar
geometry arises. At later times, if the Hill sphere becomes larger than the
scale height, polar accretion is suppressed and the flow may become
equatorial. The analytical results of this paper allow for the construction of
evolutionary scenarios of this nature, including their observational
consequences, although the latter must be studied in greater detail.

## Acknowledgements

We thank the two anonymous reviewers for their helpful comments. We thank
Andrew Bailey, Leia Barrowes, Konstantin Batygin, Leonardo Krapp, Kaitlin
Kratter, Kevin Napier, Andrew Youdin, and Zhaohuan Zhu for useful discussions
and suggestions. A.G.T. acknowledges support from the Fannie and John Hertz
Foundation and the University of Michigan’s Rackham Merit Fellowship Program.
F.C.A. is supported in part by the Leinweber Center for Theoretical Physics at
the University of Michigan. This paper made use of the Julia programming
language (Bezanson et al., 2017), the plotting package Makie (Danisch and
Krumbiegel, 2021), and Wolfram Mathematica (Wolfram Research, Inc., 2022).

## References

  * Abramowitz and Stegun (1972) Abramowitz, M., Stegun, I.A., 1972.  Handbook of Mathematical Functions.  Dover Books.
  * Adams and Batygin (2022) Adams, F.C., Batygin, K., 2022.  Analytic Approach to the Late Stages of Giant Planet Formation.  ApJ 934, 111.  doi:10.3847/1538-4357/ac7a3e, arXiv:2206.10089.
  * Adams and Shu (1985) Adams, F.C., Shu, F.H., 1985.  Infrared emission from protostars.  ApJ 296, 655--669.  doi:10.1086/163483.
  * Adams and Shu (1986) Adams, F.C., Shu, F.H., 1986.  Infrared Spectra of Rotating Protostars.  ApJ 308, 836.  doi:10.1086/164555.
  * Ayliffe and Bate (2009) Ayliffe, B.A., Bate, M.R., 2009.  Circumplanetary disc properties obtained from radiation hydrodynamical simulations of gas accretion by protoplanets.  MNRAS 397, 657--665.  doi:10.1111/j.1365-2966.2009.15002.x, arXiv:0904.4884.
  * Bate et al. (2003) Bate, M.R., Lubow, S.H., Ogilvie, G.I., Miller, K.A., 2003.  Three-dimensional calculations of high- and low-mass planets embedded in protoplanetary discs.  MNRAS 341, 213--229.  doi:10.1046/j.1365-8711.2003.06406.x, arXiv:astro-ph/0301154.
  * Batygin and Morbidelli (2020) Batygin, K., Morbidelli, A., 2020.  Formation of Giant Planet Satellites.  ApJ 894, 143.  doi:10.3847/1538-4357/ab8937, arXiv:2005.08330.
  * Benisty et al. (2021) Benisty, M., Bae, J., Facchini, S., Keppler, M., Teague, R., Isella, A., Kurtovic, N.T., Pérez, L.M., Sierra, A., Andrews, S.M., Carpenter, J., Czekala, I., Dominik, C., Henning, T., Menard, F., Pinilla, P., Zurlo, A., 2021.  A Circumplanetary Disk around PDS70c.  ApJ 916, L2.  doi:10.3847/2041-8213/ac0f83, arXiv:2108.07123.
  * Bezanson et al. (2017) Bezanson, J., Edelman, A., Karpinski, S., Shah, V.B., 2017.  Julia: A fresh approach to numerical computing.  SIAM review 59, 65--98.  URL: https://doi.org/10.1137/141000671.
  * Blandford and Payne (1982) Blandford, R.D., Payne, D.G., 1982.  Hydromagnetic flows from accretion disks and the production of radio jets.  MNRAS 199, 883--903.  doi:10.1093/mnras/199.4.883.
  * Cassen and Moosman (1981) Cassen, P., Moosman, A., 1981.  On the formation of protostellar disks.  Icarus 48, 353--376.  doi:10.1016/0019-1035(81)90051-8.
  * Chevalier (1983) Chevalier, R.A., 1983.  The enviroments of T Tauri stars.  ApJ 268, 753--765.  doi:10.1086/160997.
  * Cimerman et al. (2017) Cimerman, N.P., Kuiper, R., Ormel, C.W., 2017.  Hydrodynamics of embedded planets’ first atmospheres - III. The role of radiation transport for super-Earth planets.  MNRAS 471, 4662--4676.  doi:10.1093/mnras/stx1924, arXiv:1707.08079.
  * D’Angelo et al. (2002) D’Angelo, G., Henning, T., Kley, W., 2002.  Nested-grid calculations of disk-planet interaction.  A&A 385, 647--670.  doi:10.1051/0004-6361:20020173, arXiv:astro-ph/0112429.
  * Danisch and Krumbiegel (2021) Danisch, S., Krumbiegel, J., 2021.  Makie.jl: Flexible high-performance data visualization for julia.  Journal of Open Source Software 6, 3349.  URL: https://doi.org/10.21105/joss.03349, doi:10.21105/joss.03349.
  * Draine and Lee (1984) Draine, B.T., Lee, H.M., 1984.  Optical Properties of Interstellar Graphite and Silicate Grains.  ApJ 285, 89.  doi:10.1086/162480.
  * Fung et al. (2014) Fung, J., Shi, J.M., Chiang, E., 2014.  How Empty are Disk Gaps Opened by Giant Planets?  ApJ 782, 88.  doi:10.1088/0004-637X/782/2/88, arXiv:1310.0156.
  * Fung et al. (2019) Fung, J., Zhu, Z., Chiang, E., 2019.  Circumplanetary Disk Dynamics in the Isothermal and Adiabatic Limits.  ApJ 887, 152.  doi:10.3847/1538-4357/ab53da, arXiv:1909.09655.
  * Ghosh and Lamb (1978) Ghosh, P., Lamb, F.K., 1978.  Disk accretion by magnetic neutron stars.  ApJ 223, L83--L87.  doi:10.1086/182734.
  * Krapp et al. (2022) Krapp, L., Kratter, K.M., Youdin, A.N., 2022.  The 3D Dust and Opacity Distribution of Protoplanets in Multifluid Global Simulations.  ApJ 928, 156.  doi:10.3847/1538-4357/ac5899, arXiv:2110.02428.
  * Kuwahara and Kurokawa (2024) Kuwahara, A., Kurokawa, H., 2024.  Analytic description of the gas flow around planets embedded in protoplanetary disks.  A&A 682, A14.  doi:10.1051/0004-6361/202347530, arXiv:2311.01943.
  * Lambrechts and Lega (2017) Lambrechts, M., Lega, E., 2017.  Reduced gas accretion on super-Earths and ice giants.  A&A 606, A146.  doi:10.1051/0004-6361/201731014, arXiv:1708.00767.
  * Lambrechts et al. (2019) Lambrechts, M., Lega, E., Nelson, R.P., Crida, A., Morbidelli, A., 2019.  Quasi-static contraction during runaway gas accretion onto giant planets.  A&A 630, A82.  doi:10.1051/0004-6361/201834413, arXiv:1907.06362.
  * Lubow et al. (1999) Lubow, S.H., Seibert, M., Artymowicz, P., 1999.  Disk Accretion onto High-Mass Planets.  ApJ 526, 1001--1012.  doi:10.1086/308045, arXiv:astro-ph/9910404.
  * Machida et al. (2008) Machida, M.N., Kokubo, E., Inutsuka, S.i., Matsumoto, T., 2008.  Angular Momentum Accretion onto a Gas Giant Planet.  ApJ 685, 1220--1236.  doi:10.1086/590421, arXiv:0801.3305.
  * Marleau et al. (2022) Marleau, G.D., Aoyama, Y., Kuiper, R., Follette, K., Turner, N.J., Cugno, G., Manara, C.F., Haffert, S.Y., Kitzmann, D., Ringqvist, S.C., Wagner, K.R., van Boekel, R., Sallum, S., Janson, M., Schmidt, T.O.B., Venuti, L., Lovis, C., Mordasini, C., 2022.  Accreting protoplanets: Spectral signatures and magnitude of gas and dust extinction at H $\alpha$.  A&A 657, A38.  doi:10.1051/0004-6361/202037494, arXiv:2111.12099.
  * Marleau et al. (2023) Marleau, G.D., Kuiper, R., Béthune, W., Mordasini, C., 2023.  The Planetary Accretion Shock. III. Smoothing-free 2.5D Simulations and Calculation of H$\alpha$ Emission.  ApJ 952, 89.  doi:10.3847/1538-4357/accf12, arXiv:2305.01679.
  * Marley et al. (2007) Marley, M.S., Fortney, J.J., Hubickyj, O., Bodenheimer, P., Lissauer, J.J., 2007.  On the Luminosity of Young Jupiters.  ApJ 655, 541--549.  doi:10.1086/509759, arXiv:astro-ph/0609739.
  * Martin and Lubow (2011) Martin, R.G., Lubow, S.H., 2011.  Tidal truncation of circumplanetary discs.  MNRAS 413, 1447--1461.  doi:10.1111/j.1365-2966.2011.18228.x, arXiv:1012.4102.
  * Pollack et al. (1996) Pollack, J.B., Hubickyj, O., Bodenheimer, P., Lissauer, J.J., Podolak, M., Greenzweig, Y., 1996.  Formation of the Giant Planets by Concurrent Accretion of Solids and Gas.  Icarus 124, 62--85.  doi:10.1006/icar.1996.0190.
  * Quillen and Trilling (1998) Quillen, A.C., Trilling, D.E., 1998.  Do Proto-jovian Planets Drive Outflows?  ApJ 508, 707--713.  doi:10.1086/306421, arXiv:astro-ph/9712033.
  * Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., Sedlmayr, E., 2003.  Rosseland and Planck mean opacities for protoplanetary discs.  A&A 410, 611--621.  doi:10.1051/0004-6361:20031279, arXiv:astro-ph/0308344.
  * Shakura and Sunyaev (1973) Shakura, N.I., Sunyaev, R.A., 1973.  Black holes in binary systems. Observational appearance.  A&A 24, 337--355.
  * Szulágyi et al. (2022) Szulágyi, J., Binkert, F., Surville, C., 2022.  Meridional Circulation of Dust and Gas in the Circumstellar Disk: Delivery of Solids onto the Circumplanetary Region.  ApJ 924, 1.  doi:10.3847/1538-4357/ac32d1, arXiv:2103.12128.
  * Szulágyi et al. (2019) Szulágyi, J., Dullemond, C.P., Pohl, A., Quanz, S.P., 2019.  Observability of forming planets and their circumplanetary discs II. - SEDs and near-infrared fluxes.  MNRAS 487, 1248--1258.  doi:10.1093/mnras/stz1326, arXiv:1905.03563.
  * Szulágyi et al. (2016) Szulágyi, J., Masset, F., Lega, E., Crida, A., Morbidelli, A., Guillot, T., 2016.  Circumplanetary disc or circumplanetary envelope?  MNRAS 460, 2853--2861.  doi:10.1093/mnras/stw1160, arXiv:1605.04586.
  * Szulágyi and Mordasini (2017) Szulágyi, J., Mordasini, C., 2017.  Thermodynamics of giant planet formation: shocking hot surfaces on circumplanetary discs.  MNRAS 465, L64--L68.  doi:10.1093/mnrasl/slw212, arXiv:1609.08652.
  * Tanigawa et al. (2012) Tanigawa, T., Ohtsuki, K., Machida, M.N., 2012.  Distribution of Accreting Gas and Angular Momentum onto Circumplanetary Disks.  ApJ 747, 47.  doi:10.1088/0004-637X/747/1/47, arXiv:1112.3706.
  * Ulrich (1976) Ulrich, R.K., 1976.  An infall model for the T Tauri phenomenon.  ApJ 210, 377--391.  doi:10.1086/154840.
  * Wolfram Research, Inc. (2022) Wolfram Research, Inc., 2022.  Mathematica, Version 13.2.  URL: https://www.wolfram.com/mathematica. champaign, IL, 2022.
  * Yadav and Thorngren (2017) Yadav, R.K., Thorngren, D.P., 2017.  Estimating the Magnetic Field Strength in Hot Jupiters.  ApJ 849, L12.  doi:10.3847/2041-8213/aa93fd, arXiv:1709.05676.
  * Zhu (2015) Zhu, Z., 2015.  Accreting Circumplanetary Disks: Observational Signatures.  ApJ 799, 16.  doi:10.1088/0004-637X/799/1/16, arXiv:1408.6554.
  * Zhu et al. (2018) Zhu, Z., Andrews, S.M., Isella, A., 2018.  On the radio detectability of circumplanetary discs.  MNRAS 479, 1850--1865.  doi:10.1093/mnras/sty1503, arXiv:1708.07287.