# On the structure and mass delivery towards circumplanetary discs

Matthäus Schulik, Anders Johansen, Bertram Bitsch, Elena Lega, Michiel Lambrechts

► To cite this version:

Matthäus Schulik, Anders Johansen, Bertram Bitsch, Elena Lega, Michiel Lambrechts. On the structure and mass delivery towards circumplanetary discs. Astronomy and Astrophysics - A&A, 2020, 642, pp.A187. 10.1051/0004-6361/202037556. hal-03025362

HAL Id: hal-03025362

https://hal.science/hal-03025362

Submitted on 11 Dec 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the structure and mass delivery towards circumplanetary discs

Matthäus Schulik<fnref n="1" />, Anders Johansen<fnref n="1" />, Bertram Bitsch<fnref n="2" />, Elena Lega<fnref n="3" />, and Michiel Lambrechts<fnref n="1" />

Received 22 January 2020 / Accepted 27 March 2020

## ABSTRACT

Circumplanetary discs (CPDs) that form around young gas giants are thought to be the sites of moon formation as well as an intermediate reservoir of gas that feeds the growth of the gas giant. How the physical properties of such CPDs are affected by the planetary mass and the overall opacity is relatively poorly understood. In order to clarify this, we used the global radiation hydrodynamics code FAR-GOCA with a grid structure that allows sufficient resolution of the planetary gravitational potential for a CPD to form. We then studied the gas flows and density-temperature structures that emerge as a function of planet mass, opacity, and potential depth. Our results indicate interesting structure formation for Jupiter-mass planets at low opacities, which we subsequently analysed in detail. Using an opacity level that is 100 times lower than that of the dust of the interstellar medium, our Jupiter-mass protoplanet features an envelope that is sufficiently cold for a CPD to form, and a free-fall region separating the CPD and the circumstellar disc that emerges. Interestingly, this free-fall region appears to be the result of supersonic erosion of outer envelope material, as opposed to the static structure formation that one would expect at low opacities. Our analysis reveals that the planetary spiral arms seem to pose a significant pressure barrier that needs to be overcome through radiative cooling in order for gas to be accreted onto the CPD. The circulation inside the CPD is near-Keplerian and is modified by the presence of CPD spiral arms. The same is true when we increase the planetary potential depth, which in turn increases the planetary luminosity, quenches the formation of a free-fall region, and decreases the rotation speed of the envelope by 10%. For high opacities, we recover results from the literature, finding an almost featureless hot envelope. With this work, we demonstrate the first simulation and analysis of a complete detachment process of a protoplanet from its parent disc in a 3D radiation hydrodynamics setting.

**Key words.** planets and satellites: formation – radiative transfer – hydrodynamics – planet–disk interaction

## 1. Introduction

Circumplanetary discs (CPDs) are rotationally supported discs, consisting of gas and dust, that are thought to form around massive host protoplanets. During this formation process, gas from the parent circumstellar disc (CSD) in which the planet is formed is accreted onto the Hill sphere of the planet. Inside the Hill sphere, the gas is unable to finalise its fall onto the planet as the gas cannot easily get rid of its angular momentum with respect to the planet and therefore forms a disc (Tanigawa et al. 2012). The recent discovery of a massive protoplanet with an associated mass of possibly circumplanetary material in PDS70 (Keppler et al. 2018; Haffert et al. 2019; Isella et al. 2019) has generated renewed interest in the properties of this class of objects.

Circumplanetary discs have long been recognised as potential sites for the formation of regular moon populations around giant planets (Canup & Ward 2006; Shibaike et al. 2019; Ronnet & Johansen 2020). Because they are directly coupled to the global CSD gas flows, CPDs are thought to be possible bottle-necks for further growth of giant planets. Those global flows can exhibit supersonic spiral arm shocks as the planet orbits relative to the CPD medium (Goldreich & Tremaine 1978), and once the planet grows massive enough, those shocks perturb the global flows sufficiently to open gaps in the gas distribution (Goldreich & Tremaine 1980). These gaps are thought to further limit the accretion rates onto the planetary Hill spheres and in turn are

thought to regulate CPD and planet growth (D'Angelo & Lubow 2008).

On the other hand, radiative feedback from the growing planet and CPD into the CSD is inevitable, with potentially observable consequences for the CSD chemistry (Cleeves et al. 2015). The study of such a triplet system of planet-CPD-CSD must therefore be performed numerically in order to access all available physical information.

Because of the inherent complexity of simulating this problem, hydrodynamical simulations are usually either limited in the number of CSD orbits they can compute, or in the physics that is being computed. In this endeavour, a number of interesting milestones exist: D'Angelo et al. (2003) used the system of 2D Euler equations with a cooling prescription and found various CPD structures with spirals in them. Contrary to this latter study, Ayliffe & Bate (2009a) found no CPD spiral arms in a global 3D radiation hydrodynamics setting, arguing that they disappear due to the increase in degrees of freedom for the flow in 3D. Tanigawa et al. (2012) found a disc in Keplerian rotation in their isothermal 3D simulation. Gressel et al. (2013) showed in a comparison of global isothermal, adiabatic-hydrodynamic and MHD simulations that isothermal simulations favour highly Keplerian disc rotation and strong flattening of the CPDs. Zhu et al. (2016) again found spiral arms in highly resolved radial-azimuthal 2D runs with a cooling prescription as a function of optical depth, identified them to be triggered by the tidal action from the host

<footnote n="1">¹ Lund Observatory, Box 43, Sölvegatan 27, 22100 Lund, Sweden
e-mail: schulik@astro.lu.se</footnote>

<footnote n="2">² Max-Planck Institut für Astronomie, Königsstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="3">³ Laboratoire Lagrange, UMR7293, Université de la Côte d'Azur, Boulevard de la Observatoire, 06304 Nice Cedex 4, France</footnote>
<table><caption>**Table 1.** Simulation runs used in this paper.</caption><thead><tr><th>Label</th><th><i>m</i><sub>P</sub>/<i>m</i><sub>Jup</sub></th><th><i>m</i><sub>P</sub>/<i>m</i><sub>⊕</sub></th><th><i>r</i><sub>s</sub></th><th>Opacities <sup>(a)</sup></th></tr></thead><tbody><tr><td colspan="5">CPD occurrence</td></tr><tr><td>m1</td><td>0.06</td><td>20</td><td>0.1</td><td>1.0/0.01</td></tr><tr><td>m2</td><td>0.2</td><td>60</td><td>0.1</td><td>0.01</td></tr><tr><td>m3</td><td>0.3</td><td>90</td><td>0.1</td><td>1.0/0.01</td></tr><tr><td>m4</td><td>0.4</td><td>120</td><td>0.1</td><td>1.0/0.01</td></tr><tr><td>m5</td><td>0.75</td><td>240</td><td>0.1</td><td>1.0/0.01</td></tr><tr><td colspan="5">Main simulations</td></tr><tr><td>C1 (“nominal”)</td><td>1</td><td>360</td><td>0.1</td><td>0.01</td></tr><tr><td>H1 (“high opacity”)</td><td>1</td><td>360</td><td>0.1</td><td>1.0</td></tr><tr><td>H2 (“deep”)</td><td>1</td><td>360</td><td>0.025</td><td>0.01</td></tr><tr><td>C2 (“Bell & Lin”)</td><td>1</td><td>360</td><td>0.1</td><td>0.01</td></tr></tbody></table>

**Notes.** The relation of the runs with the science questions posed in this work is explained in Sect. 1 and parameters are explained in Sect. 2. <sup>(a)</sup>Opacities used: numbers denote constant opacities $κ$ for all simulations in cm$^2$ g$^{-1}$, except for simulation C2, which uses Bell & Lin opacities reduced by the factor in the table. All simulation runs have an unperturbed CSD surface density of $\Sigma_0 = 100$ g cm$^{-2}$ at the position of the planet. Simulation letters refer to the purpose or expected simulation outcome: m = “mass survey”, C = “cold”, H = “hot”.

star, and measured accretion rates onto the planet through the CPD. The same study also investigated the effects of the optical depth on the CPD properties in detail by changing the CPD mass. An interesting result from the approach of these latter authors of varying the CPD mass is that the temperature changes with the optical depth. The temperature in turn controls the attack angle of the spiral arms, which in their study is responsible for the ongoing accretion of the CPD onto the planet. Subsequently, Szulágyi et al. (2016) clarified in 3D radiation hydrodynamics runs that the occurrence of CPDs, as opposed to spherical envelopes around giant planets, is linked to their temperature; that is, under otherwise identical conditions a hot envelope will collapse into a disc when cooled artificially. In total, those studies made clear that the exact accretion rates and the rotational and structural properties of simulated CPDs depend on the numerical framework and the treatment of thermodynamics.

Additionally, Schulik et al. (2019, hereafter S19) showed that underresolving the central planetary region raises temperatures (and entropies) to incorrect levels. Furthermore, reducing the opacities of the planetary envelope was demonstrated by these latter authors to significantly aid the flattening of envelopes surrounding Saturn-mass planets. Therefore, there is a clear need to re-visit past results concerning CPDs in the framework of 3D global radiation hydrodynamics simulations with sufficient numerical resolution and a variation of opacities for Jupiter-mass planets.

With this motivation, we ask and address a number of science questions related to our two sets of simulations listed in Table 1:

– At which protoplanet masses do CPDs appear? We address this by quantifying the flatness and state of rotation of circum-planetary envelopes for a number of protoplanets with increasing masses in runs m1-m5, and C1 / H1.

– What are the masses and accretion rates of CPDs relative to those of the planet? Due to the coupling of gas flows from the CSD into the planetary Hill sphere, CPDs will to some degree co-evolve with their host planet. This process will be particularly important if the CPD is accreting onto the planet and at the same time is replenished from the CSD. Earlier authors have

often assumed that the mass accreted by the planet is a fraction of the mass entering either the Bondi or the Hill radius (D'Angelo & Lubow 2008; Owen & Menou 2016), but a more recent study (Lambrechts et al. 2019) shows that the mass flux through the Hill sphere and the actual accretion rate can differ by a factor of up to 100, with this number approaching ~0.1 when the planetary mass reaches ~1 $m_{Jup}$. We use runs C1/C2 and H1/H2 to further investigate this phenomenon and measure the mass and accretion rates of the CPDs and compare them to those of the planets.

– What is the influence of the variation of dust opacity on properties of the CPD? Is there an influence from the CSD on the accretion process? Discs at different ages are able to provide differing amounts of dust to already formed protoplanetary envelopes (Brauer et al. 2008; Birnstiel et al. 2012). However, the full-time evolution of CPDs in 3D radiative hydrodynamics simulations is currently not accessible on modern hardware. This problem is even more severe if one would consider the evolution of the full disc, including stellar irradiation and the evolution of the size distribution of the dust. We therefore take a simplified approach and compare a dust-rich (run H1) to a dust-poor (run C1) scenario. We also perform two support parameter simulations, investigating deep planetary potentials (run H2) and more realistic, non-constant opacities provided by Bell & Lin (1994) (run C2).

We structure this paper as follows in order to address the above questions. In Sect. 2, we describe our methods, the used set of simulations and justify parameter choices. Section 3 describes how the envelopes of protoplanets of varying masses and envelope opacities flatten and rotate. Section 4 describes and analyses our nominal simulation run, which presents most of the novel and interesting physics features. Subsequently, Sect. 5 discusses how those features change when using different simulation parameters and establishes the robustness of our nominal results. Section 5 also presents results for CPD masses, mass profiles, and accretion rates, together with a related discussion.

## 2. Methods

The methods used to address the science goals are mostly identical to our previous work, namely Schulik et al. (2019, S19). We briefly reiterate these methods here and point out major changes.

### 2.1. Physical model and parameters

We use the code FARGOCA (Bitsch et al. 2014; Lega et al. 2014) to solve the equations of viscous hydrodynamics coupled with an evolution equation for the zeroth moment of the photon intensity, which is the mean photon energy density equivalent to a separate photon temperature. The hydrodynamic advection problem is solved using the FARGO algorithm (Masset 2000), while momentum transport evolves also via the full viscous tensor through the employment of a classical Strang-splitting scheme. The radiation transport is solved with the help of flux-limited diffusion (Levermore & Pomraning 1981) and we set $κ_R = κ_P ≡ κ$ for reasons of simplicity. The equation of state for the gas is adiabatic with a constant adiabatic index of $γ = 1.4$. The mean molecular weight is that of a solar hydrogen-helium mixture, which is $μ = 2.35$ g mole$^{-1}$. Viscous momentum transport is computed through the full hydrodynamic stress tensor multiplied with a constant physical viscosity of $ν = 10^{-15}$ cm$^2$ s$^{-1}$ that corresponds roughly to $α ≈ 10^{-2}$ at our disc temperatures. Fujii et al. (2011, 2014) found that CPDs will not be hot enough to be ionised, and therefore we do not include the effects of magnetohydrodynamics in our simulations. The main parameters of
interest for our work are the planetary masses, the opacities, and the gravitational smoothing lengths. Their variations and corresponding simulation labels are listed in Table 1.

The planetary masses used are $20 m_⊕$ up to $360 m_⊕$; the latter is just slightly above $\approx 1 m_{Jup}$. We use two different constant opacities, as listed in Table 1, to scan the parameter space of CPD properties. High-opacity runs ($\kappa = 1.0 \text{ cm}^2 \text{ g}^{-1}$) serve to connect our study to previous work by Lambrechts & Lega (2017) (for the $20 m_⊕$ case) as well as Szulágyi & Mordasini (2017) and Lambrechts et al. (2019) (for the Jupiter-mass planets).

Quantities denoted with a tilde, such as $\tilde{x} = x/x_H$, have numerical values indicating the Hill radius fraction. This is similar to the notation used in Tanigawa et al. (2012), where all quantities were normalised with the pressure scale height $h$, but as these latter authors used $r_H = h$, the locations of circulation features in our plots can be directly compared to those in their study.

Gravity in our simulation domain is determined by contributions from both the star and the planet. The planetary gravity consists only of a gravitational potential with an inner cut-off defined by the smoothing length $r_s$ (Klahr & Kley 2006). This is a key numerical parameter in our work. All our simulations are performed with well-resolved smoothing lengths (ten cells per $r_s$, as defined in S19) and a nominal potential depth of $\tilde{r}_s = 0.1$. This value was chosen after reviewing the literature on CPDs, indicating CPD outer edges at around $\tilde{r} = 0.3-0.5$. In one case, that of our most expensive simulation (run H2), we go as deep as $\tilde{r}_s = 0.025$, which corresponds to a physical size of $r = 12 r_{Jup}$ at our orbital distance of 5.2 AU.

The independent evolution of gas internal energy and mean photon energy density allows the code to find regions where those two quantities are either at equilibrium or are uncoupled from each other. The former is the case in dense, optically thick media while the latter is the case for very tenuous regions of the simulation domain. The coupling between gas and photon energy is calculated using the Planck-mean opacity, while the diffusion and/or escape of photons into space is calculated using the Rosseland-mean opacity; both are taken as equal in this work and we subsequently only refer to them as the opacity. We then use either different values of constant opacities $\kappa$ or those from Bell & Lin (1994) scaled by a factor, which we list as $\epsilon$ in the last column in Table 1. Those opacity values are compatible with the upper and lower limits of opacities suggested in Mordasini (2014) which correspond to a typical factor of $10^2$ to $10^4$ reduction compared to interstellar medium (ISM) dust opacities. The exact equations that we solve for the radiative transfer can be found in S19. Planetary gaps are generated via the same process as in S19 and are evolved for 400 orbits before the resolution is increased to resolve the CPD.

## 2.2. Numerical setup and grid refinement

We work in a global, spherical coordinate system that is centred on the star and label the independent coordinate as $[r, \theta, \phi]$. We also use the planetocentric coordinates $[x, y, z] = [(r - 1) \cos \theta \cos \phi, r \sin \theta \cos \phi, r \sin \phi]$, so that the planet will be at $[x_p, y_p, z_p] = [0, 0, 0]$ and we often refer only to the planetocentric coordinates for simplicity.

Boundary condition data are generated as in S19 and as in Kley et al. (2009), Lega et al. (2014), Bitsch et al. (2014) via the use of three distinct numerical steps. The first generates a radiative, radial-vertical disc equilibrium in 2D with a sufficiently wide radial extent. In the second step, the gap-formation step, the data from the first step are run in low, equidistant resolution

in 3D for 400 orbits, which has been previously shown in S19 to be deep enough for significant gap depths. Finally, the third step takes a limited radial extent of the global disc and runs it in high resolution around the planet, as explained below.

The simulation domain in all our simulation runs is $\theta \in [-\pi, \pi]$, $\phi \in [81, 90]^{\circ}$ for azimuth and colatitude. Boundary conditions in azimuth are periodic for all variables. In the colatitude, we use a half disc, that is, reflective boundaries at $\phi = 90^{\circ}$ for all variables. Reflective boundaries are used for hydrodynamic variables at $\phi = 81^{\circ}$, and open boundaries for the radiative energy at $\phi = 81^{\circ}$. The radial extent of our simulation domains is adjusted depending on the planet mass. For $20 m_⊕ < m_p < 120 m_⊕$ it is sufficient to run with $r \in [0.7, 1.3] \times 5.2$ AU. However, once planetary gaps become deeper for higher masses, the gap width also increases. Hence, we run the massive protoplanets with $m_p > 120 m_⊕$ with $r \in [0.4, 1.6] \times 5.2$ AU. Using an overly narrow radial simulation for the high-mass planets results in artificial gap edge instabilities that feed large surges of gas onto the planet. These instabilities disappear once running with the larger radial extent.

The spacing of the simulation grid in radial, azimuthal, and colatitudinal dimensions is of the form $d\Xi \propto 1 - g(\Xi - \Xi_0)$, where $g(\Xi)$ is a Gaussian function, $\Xi$ represents any of the three spatial directions $r, \theta, \phi$ in our spherical coordinate system, and $\Xi_0$ is the planetary position for any of the three coordinates. This type of grid results in a highly resolved, quasi-cartesian grid in the direct vicinity of the planet. We design the three grids in such a way that the increase in resolution element size has the same value for all three spatial dimensions in order to not introduce artificial effects of over-resolution in any one direction. In order to minimise computation time, the resolution gradients are stretched to a reasonable maximum that leaves the horseshoe orbits unchanged and allows radial density gradients in the gap to converge numerically.

The choice of the minimum resolution element $d\Xi_{\min}$, that is, the element with the best resolution directly adjacent to the planet and maximum of $g(\Xi)$, is given by the requirement for the gravitational smoothing length to be resolved by ten cells and is therefore determined after the choice of $r_s$. This requirement originates from our previous studies focused on accretion rates (S19), but we find that adequately resolving $r_s$ is equally important for the structure of CPDs in order to avoid artificial overheating and artificial flattening of the planetary envelopes. For the nominal simulation run this results in a minimum resolution element of $dx_{\min} = 7 \times 10^{-4}$ (or $7.5 r_{Jup}$) in size and our simulation run H2 with the deep potential is of $dx_{\min} = 1.7 \times 10^{-4}$ (or $1.9 r_{Jup}$).

# 3. Results: occurrence of CPDs with increasing planetary mass and the influence of opacity

We now present the results from our mass survey. In general, we refer to planetary envelopes when describing the gas surrounding the smoothing length, as this nomenclature remains agnostic towards the existence of a CPD.

## 3.1. Measuring the flatness of envelopes

In order to assess the structure of the planetary envelopes, one could pursue the idea of measuring the envelope aspect ratio $H/r$. However, $H/r$ does not provide information about the Keplerian rotation support of the envelope. This is because without an assumption about the ratio of pressure to centrifugal
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Rotation profiles in the envelopes of planets of varying masses. We define the regions in the surroundings of the planet according to the rotational state of the envelope. We note that the region labelled CPD is only appropriately named for masses of $m_p \ge 240 m_\oplus$. The strict transition into the Keplerian shear of the CSD happens at $v_\theta/v_{\text{Kepler}}=-1$, as can be seen from the slope turnover in the rotation curves. The vertical black dashed curve denotes $r=0.3r_H$, where rotation and density asymmetry values are measured and plotted in Fig. 2.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Occurrence of CPDs as measured by the flatness parameter and Keplerianity of planetary envelopes of different masses. Simulation data are taken in steady state after a runtime of five orbits. The maximum value of $v_\theta/v_K$ is also indicated in order to monitor fluctuations along one orbit. While the Keplerian rotation fraction of the protoplanetary envelopes is rising quasi-linearly with mass, the flatness shows a sharp decrease at 120 $m_\oplus$ for both values of the constant opacities. This shows that the infall and accretion of angular momentum in the midplane follow a straightforward relation with the potential that they are accreted into. On the other hand, the vertical cooling responsible for the envelope flattening plays an independent role and prevents CPD formation at high opacity, even if the necessary potential depth is given.</figcaption></figure>

support, the aspect ratio can only be recast into the form

$$ \frac{H}{r} = \left( \frac{c_s}{v_k} \right)^2, \qquad (1) $$

which is not $v_\theta/v_k$. Hence, in order to measure the flatness of a disc we define the flatness parameter at distance $d$ through the cylindrically averaged 2D density $\bar{\rho}(r, z)$:

$$ \xi(d) = \frac{\bar{\rho}(0, d)}{\bar{\rho}(d, 0)}. \qquad (2) $$

This definition is a simple measure of the asymmetry between the density profiles in the vertical versus midplane directions and is thus an indicator for “disciness”. We use this parameter $\xi$ additionally to $v_\theta/v_k$ to assess the properties of planetary envelopes with increasing mass for the simulation sets m1–m5, C1, and H1.

## 3.2. Envelope flatness

We first present the simulation data which we used to inform our later decisions on identifying key interesting simulation parameters to be studied in detail. We investigated the rotational state of the midplane around planets of various masses at a time of five orbits after the start of our simulation runs. This allows the planetary envelope enough time to settle into rotational equilibrium, to the degree that the counteracting pressure gradients permit this. This is also enough time to establish an accretion equilibrium between gas flowing from the envelope onto the planet and accreted high-angular momentum gas replenishing gas in the envelope.

We define the distinction between planet, circumplanetary disc, and envelope as follows. The planet is all the material

inside the sphere of radius defined by the smoothing length, that is, $r=r_s$, centred on the centroid coordinates of the planetary potential, namely $(r_p, \theta_p, \phi_p) = (1, 0, \pi/2)$. We emphasise that this is neither a real planetary surface nor a planetary interior. We define the CPD as the region inside of $0.5r_H$ with prograde, or positive, rotation. This definition for the CPD excludes the region from $0.5 < \tilde{r} < 1.0$, as this is a region into which horseshoe orbits penetrate (seen in Fig. 1 as regions of retrograde, or negative, rotation), which we refer to as the outer envelope. The sum of CPD and outer envelope we simply refer to as the envelope. Hence, material entering the Hill sphere will enter the envelope through the outer envelope.

For all planet masses, we measure the rotation profiles in their envelopes and plot them in Fig. 1 normalised to the individual Hill radius of each planet. Between planetary masses of 20 $m_\oplus$ and 60 $m_\oplus$ we see an important evolution of the rotational profile. A general steepening and retreat of the transition into the Keplerian shear of the disc is evident. From 60 $m_\oplus$ to 120 $m_\oplus$ there is only a relatively weak evolution in the envelopes. After a mass-doubling from 120 $m_\oplus$ towards 240 $m_\oplus$, a transition in the shape of the rotational profile occurs, developing a “pedestial” consisting of relatively high Keplerian rotation in the CPD region for the 240 $m_\oplus$ and 360 $m_\oplus$ planets.

Since the parameter space of the envelopes of planets around the classical, critical runaway core mass have already been studied in detail in 3D radiation hydrodynamical settings (Ormel et al. 2015; Lambrechts & Lega 2017; Kurokawa & Tanigawa 2018), here we focus on the high-mass end of planet evolution. A key parameter regulating the ability of a gaseous envelope to rotate is the opacity of the gas–dust mixture. In order to assess its importance in forming rotating CPDs, our scan in planet masses was performed with two different constant opacities.
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Overview of flow structures in the midplane of the nominal simulation run C1 after reaching steady state at orbit 5. Density (**left**) and temperature (**right**) values in the midplane are shown along with the gas streamlines coloured according to Mach number. While streamlines in all generality do not coincide with the gas motion, they do so in steady-state. Gas coming from the circumstellar disc encounters the spiral arms in the midplane, which are heavily modified by the ongoing accretion process and non-isothermal modification of the planetary Hill sphere. The rotationally supported disc is formed mainly between $\tilde{r}=0.1$ and $\tilde{r}=0.4$. This disc contains its own structures, such as a midplane accretion shock from the supersonically infalling gas and smaller CPD spiral arm shocks. Arrows target the features, not specific points. The letter-labels refer to the CSD spiral arm (a), free-fall region (b), CPD accretion shock (c), and CPD spiral arm (d). The black stars denote the positions where streamlines of interest cut through the spiral arm shock. These positions are constant in 3D-space, and will serve as orientation points when investigating the vertical direction. We note the highly supersonic flows which remain supersonic after encountering the spiral arm shock. This is due to insufficient static shock pressure, which is advecting the spiral arm radially inwards. The spiral arm mass is replenished from the vertical direction.

The flatness at $r = 0.3 r_H$, being representative of the strongly rotating part of the envelope, and the fraction of Keplerian rotation after five orbits in steady state are plotted in Fig. 2. It is evident that at higher masses in particular, opacity plays an important role in setting the flatness and the rotational state of envelopes. The curve for the maximal values of $v_\theta/v_{\text{Kepler}}$ indicates the deviation of individual fluid elements along the orbit of radius $0.3 r_H$ and shows that the CPD for the Jupiter-mass planet with constant opacity of $\kappa = 0.01 \text{ cm}^2 \text{ g}^{-1}$ rotates non-uniformly with individual fluid elements reaching up to 90% of the Keplerian value.

Figure 1 clearly shows that protoplanets must reach approximately Jupiter-mass before a significant rotationally supported CPD forms. We therefore focus in the following section on analysing the gas flow and density–temperature structure around the Jupiter-mass planet in the nominal simulation C1, which exhibits a richness of physical features. In Sect. 5 we continue to compare the C1 structure to the other Jupiter-mass simulations where we vary the physical and numerical parameters (runs C2-H2).

**4. Results: gas flows and envelope structure in the nominal simulation run**

This section presents the features that we find in the nominal simulation run C1. We first show and explain the simulation outcome in the midplane. We then do the same for the vertical direction, which exhibits more complex physics that couple into the mid-plane dynamics. For both midplane and the vertical direction, we describe three separate topics: the density and temperature structure of the envelope, mass delivery towards the CPD, and circulation inside the CPD.

**4.1. The midplane structure**

In Fig. 3, we show the density and temperature structure together with the associated midplane flows. The innermost region

around the planet forms an approximately spherically symmetric density and temperature distribution that extends beyond the smoothing length at $\tilde{r} = 0.1$ towards $\tilde{r} = 0.2$. A radial inflow exists which slows down as the pressure support increases. Between $\tilde{r} = 0.1$ and $\tilde{r} = 0.2$ the pressure support of the envelope gradually fades until at $\tilde{r} = 0.2$ the centrifugal forces become sufficient to flatten out the envelope and let the gas rotate with 80–90% $v_K$; see Fig. 1.

Between $\tilde{r} = 0.2$ and $\tilde{r} = 0.4$ we find the CPD proper, with time-independent spiral arm features that are evident in density as discontinuity, and in temperature as local temperature maximum (marked as feature (d)). Here, we refer to these spiral arms as “CPD spirals”. The CPD spirals share a superficially similar morphology with those presented in Zhu et al. (2016). However due to our limitations in resolution and gravitational smoothing, we only see the outermost region of what Z16 are able to probe. Compared to this latter study, our spiral arms are relatively thick, which we attribute to our high viscosity, while the simulations of this latter study are inviscid. We discuss the three-dimensional effects of the spiral arms in Sect. 4.8, where we discuss 3D flows inside the CPD.

At $\tilde{r} \approx 0.4$ we find a midplane accretion shock, which results from the collision of the free-falling gas and the CPD (feature (c)). Its effects can be seen as a sharp outer rim in the midplane temperature at the CPD edge. Between $\tilde{r} = 0.4$ and $\tilde{r} = 0.9$ there is a region of free-falling CSD gas, which exhibits an inversion in pressure gradient (feature (b)). We subsequently name this region simply the free-fall region (FFR). Correspondingly, this region is evacuated relative to the CSD region and the density and temperature decrease inwards in the midplane as the radial velocity increases, which we plot in Fig. 4.

At even larger distances than $\tilde{r} \approx 0.9$ the spiral arms shock (feature (a) in Fig. 3), decelerate, and funnel material into the FFR, similarly to Bondi-Hoyle type accretion, as opposed to Kelvin-Helmholtz accretion which is regulated by cooling and contraction of the envelope. This interaction between the CSD
<figure><img src="image_5.png" /><figcaption>Fig. 4. Overview of slices through the structures for the non-circled star in the nominal run C1 at orbit 5 for the midplane. The letter labels correspond to the same features as in Fig. 3, only on the opposite side. In the midplane, the most prominent feature is the free-fall region (marked with b and braces, grey area) just after the spiral arm shock (a). The free fall is notable as the density profile decreases as the radial velocity accelerates. The free fall is terminated when the gas hits the CPD at the left edge of the grey area, seen as a sharp increase in density and as CPD accretion shock (c). We note that we use a linear density scale here for emphasis. The CPD spiral arm can be seen as a temperature bump (d). The black star denotes the same position as in Fig. 3.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Cuts through the pressure support vs. radius inside the Hill sphere for the same times as Fig. 6. Of particular interest is the state of the simulation at $t=0 \Omega$ just after the gap formation run. We find the envelope is in a state of latent imbalance due to the low resolution in the gap formation run. From there, the CSD flows push into the envelope until a new equilibrium is found. The inset shows the evolution of the static shock pressure support. The difference in position between the star and the shock positions at $t=0-2\Omega^{-1}$ showcase the slow evolution of the spiral arms between those times. The black star denotes the same position as in Fig. 3.</figcaption></figure>

spiral arm shocks, the accretion flow, and the physical origin of the FFR are further discussed in Sect. 4.4.

Moving to the largest radial distances of interest, in the range $0.9 < \tilde{r} \lesssim 6$, the planetary gap in the CSD becomes optically thin. This implies that photons emitted from the planet do not heat its immediate vicinity, but are only re-absorbed at the gap edge, keeping the gas flowing into the Hill sphere cool at $\approx 30$ K.

## 4.2. Free-fall onto the CPD in the midplane

We now turn to a more detailed investigation into the origin of the free-fall region (feature (b)). A profile of the state of the hydrodynamic variables after the opening of the free-fall region can be seen in Fig. 4. In order to investigate the opening process of the FFR, we use the data from three snapshots in time taken at $t=0, 1, 2 \Omega^{-1}$. For those times, we plot the pressure gradient normalised to the local gravity in the same slice through the non-circled star from Fig. 3, which we present in Fig. 5.

There, it is evident that the initial pressure support in the outer envelope ($\tilde{r} \approx 1$) is near zero. This allows the circumstellar flows to penetrate the envelope (seen in the inset, as the shock pressure evolves upwards). A region of negative pressure support, i.e. free-falling gas, is subsequently established in the outer envelope.

We expand this analysis by plotting the evolution of the density together with the streamline velocities in the midplane in Fig. 6. In this figure, it becomes evident that the opening process seems to commence at a distance of around $1r_H$ just behind the spiral arm, and later continues to empty out an entire region behind the spiral arm in an asymmetric fashion. We think that efficient advection is responsible for this process. Advection is naturally more efficient for faster gas flow, and therefore the erosion of gas that opens the free-fall region is facilitated when the gas remains supersonic after the shock. Remaining supersonic after the shock is only possible for the flow when the ram

pressure of the incoming CSD material overwhelms the shock ram pressure. Therefore, in order to analyse this process further, we define the ratio of total pressure pre-shock to static pressure post-shock as

$$ \zeta \equiv \frac{(p_{\text{dyn}} + p_{\text{stat}})_{\text{pre}}}{(p_{\text{stat}})_{\text{post}}} = \frac{(\rho u_{\perp}^2 + \rho \gamma c_s^2)_{\text{pre}}}{(\rho \gamma c_s^2)_{\text{post}}}, \quad (3) $$

where $u_{\perp}$ corresponds to the streamline-velocity component perpendicular to the shock surface. We stress at this point that the nominator and denominator of $\zeta$ are generally independent of each other due to the supersonic nature of the inflow from the circumstellar disc. We use the definition of $\zeta$ to analyse the time evolution that is plotted in Fig. 6. There, streamlines that are decelerated to subsonic post-shock speeds can be identified as having $\zeta < 1$. Those that remain supersonic post-shock coincide with the area at the shock where $\zeta > 1$. We continue using this quantity to greater extent in the vertical direction, but are already able to draw some conclusions: once the incoming CSD flows possess sufficient dynamic pressure to overwhelm the static shock pressure, those flows react to the shock as if it were a “speed bump” but not a dominant counteracting force, and continue at supersonic speeds. This helps enormously in transforming the initially massive spiral arm into a low-density region between two split, weaker arms. If the static shock pressure is high enough to act against the incoming flow, then the corresponding streamlines behave like a text-book shock, being decelerated according to the jump allowed by the Rankine-Hugoniot conditions in the planetary co-rotating rest frame. This can also be viewed from a different perspective. At sufficiently high Mach-number, the post-shock material can always be advected with the flow. This would however impart a nonzero shock velocity. In the rest-frame of this shock, the post-shock material will always be subsonic, which we confirmed using a toy shock model. The remaining puzzling issue therefore
<figure><img src="image_7.png" /></figure>

**Fig. 6.** Evolution of the FFR as seen in velocity as streamline colour and density as background colour, with density contours to guide the eye. Snapshots are taken at intervals of 1.0 Ω⁻¹. Symbols along the spiral arm shock surface denote the pressure ratio ζ, as defined in Eq. (3). Filled circles denote approximate ram-pressure equilibrium, i.e. ζ ≈ 1, and filled triangles indicate ζ > 1, i.e. the region where the streamlines can push past the shock. The ζ > 1 streamlines remain supersonic after encountering the shock. Those post-shock supersonic streamlines evacuate the envelope efficiently, which helps accelerate the flow further (note the intensification of the red colour on the streamlines between the *middle* and *right panels*, which also corresponds to an increase of velocity in absolute numbers). We note also how the original spiral arm splits in two through this process.

pertains to how the spiral arm can remain in existence, as seen in Fig. 6. We show in the following sections that the vertical mass flux is most likely responsible for this.

It is important to note that this process has to take place close to the planet. Otherwise, one can always find streamlines that satisfy ζ > 1, as through the Keplerian shear in the CSD $ρu^2$ increases more than linearly. Furthermore, most of the streamlines far away from the planet at the gap edge fulfil this criterion.

It therefore seems that the opening of a free-fall region is permitted by the closeness of the ζ > 1 streamlines to the planet. This closeness to the planet is a result of the temperature, which determines the transition ζ = 1. Prior to the opening process, temperatures in the soon-to-be free-fall region are as cold as 75 K, which is even colder than assumed in some isothermal simulations, for example those of Zhu et al. (2016) where the midplane temperature inside $r_H$ is truncated to 100 K.

A last comment is in order to clarify the nature of the free-fall region in relation to numerical solutions found in other works. One might expect a planet to form a relative vacuum between the disc and itself, as for example in Béthune (2019) for strongly contracting planets, or as static solution for a given mass at low constant opacity, as ours. However, as the above analysis reveals, the mechanism for opening a free-fall region appears to be a dynamic process in conjunction with the thermodynamics of the spiral arms, rather than the properties of low-dimensional static solutions.

### 4.3. Rotation of the CPD and CPD spiral arms in the midplane

The flow that enters the CPD through the midplane triggers the midplane accretion shock at around 50% of all midplane angles. However, with only 50 K above the 200 K CPD background temperature, this shock is fairly weak and presumably has no significant influence on the entropy evolution of the CPD. After the accretion shock, the gas quickly encounters the CPD spiral arms. Those spiral arms torque the flow significantly, but in such

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

**Fig. 7.** Midplane rotational structure for the run C1 with $|v|/v_{Kep}$ as contours. After gas from the FFR enters the CPD, it encounters one of the spiral arms and is shocked towards higher values of $|v|/v_{Kep}$. This results in a de-facto isolation between on one hand the infalling material into the FFR and on the other hand the planet in the midplane. This material nevertheless does find its way into the planet through a complex vertical circulation, as is clear from the streamlines in Fig. 10. The blue, dashed circle marks the gravitational smoothing length.

way as to transfer radial momentum to angular momentum, and thus boost the flow to higher fractions of $v_θ/v_K$ than those it initially possesses. This can be observed in Fig. 7 where, following the streamlines, a region of high Keplerian rotation is evident after the encounter with the CPD spiral arms. This behaviour of the spiral arms is distinctly different from that seen in Zhu et al. (2016), where the CPD spiral arms act to reduce the angular momentum of the flow.

The influence of spiral arms on the mean flow can be understood qualitatively in terms of interpreting the spirals as oblique shocks. An oblique shock is a shock that has an inclination with respect to the local flow. These have already been studied for a
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Overview of density structures in a vertical cut through the non-circled star-symbol along the *y* = 0-axis in Fig. 3 for the nominal simulation run C1. The red circle marks the approximate Hill sphere. An approximate boundary of the CPD is the vertically integrated optical thin–thick transition, marked as τ = 1. Features are labelled identically to the temperature plot (below) with letters as follows: classical spiral arm (a), midplane free-fall region (b), (midplane) accretion shock into the CPD (c), CPD spiral arm (d), vertical accretion shock (e), and accretion funnel from colliding flows (f).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Overview of temperature structures for the same plane as in Fig. 8. Some features are more clearly distinguishable than in the density plot (above) in particular: the midplane accretion shock (c) and the vertical extension of the CPD spiral arms (d). Those features exhibit a radial asymmetry because the direction of the cut along the *y* = 0-axis does not coincide with the symmetry axis of the CPD.</figcaption></figure>

long time in the hydrodynamics literature (Kevlahan 1997, and references therein) where it is made clear that there is a vorticity jump imposed on the local flow through the shock. The vorticity jump is in general given by the density jump and the shock curvature with respect to the local flow. We note that because the CPD spiral arm shock is very “fluffy” we cannot analyse the shock in an equally quantitative manner to that with which we analyse the CSD spiral arms in the following section, and therefore a fit of the post-shock values fails for the CPD spirals.

The inner CPD boundary also merits further discussion. Here, gas does not flow directly into the planet, although it orbits at sub-Keplerian velocity. This is partially an effect of substantial pressure support. On the other hand, gas that originates in the CPD still flows into the planet, although in a vertical manner, which is explored further in Sect. 4.8.

4.4. Vertical structure

The vertical structure in density and temperature is plotted in Figs. 8 and 9, which are discussed in this section. There, vertical extensions of phenomena in the midplane can be identified. We additionally plot the vertically integrated optically thin–thick transition marked with “τ = 1”, as it approximately delineates the boundary of the CPD.

Feature labels in those figures are identical to those used pre-
viously in Figs. 3 and 4, and are additionally described with (r, z)

coordinates in the text for clearer identification. Figures 8 and 9
have a radial extent chosen such that the density plot highlights
the gap density gradient, and their vertical extent is chosen such
that the tilted structure of the spiral arms at our chosen lower
density cutoff is still visible.

The CSD spiral arms are tilted in the vertical direction (feature (a)). As we show below, this is a consequence of the disc thermodynamics-related process discussed earlier for the FFR. In the midplane, directly inward from the spiral arms, the FFR gas falls (feature (b)) onto the CPD. The free-falling gas from the CSD in the midplane can be best identified in the density plot, where it emanates from (x, z) = (±0.9, 0) at the position of the spiral arm shocks. From those positions it spreads out into a vertical fan until it hits the CPD between (x, z) = (±0.4, 0) and (x, z) = (±0.6, 0.2). At those positions, the vertical extent of the CPD accretion shock is visible as a temperature maximum that traces a concave path (feature (c)).

Next, we find the CPD spiral arms and their vertical extent (feature (d)). The CPD spiral arms show the best contrast in the temperature data, and emanate from (x, z) = (±0.3, 0) towards higher altitudes and then towards the vertical accretion shock in a looped structure. We note that this loop-structure (particularly clearly seen in Fig. 10) does not result from any vertical compression, as the gas flows parallel with respect to the loop. We also observe that this structure shows a time-variance that is identical to the tidal arms in the midplane. From this, we propose that the
<figure><img src="image_11.png" /></figure>

**Fig. 10.** Overview of a slice of the gas velocity field in the same plane as in Figs. 8 and 9. Arrows correspond to the radial and vertical components of the velocity on the slice and the colour shows the norm of the full 3D velocity. Most of the gas accreted by the planet enters the Hill sphere through the midplane and the tilted spiral arm shocks. The spiral arms give vertical kicks to the passing flows, forcing them to flow along the shock downwards into the midplane, where a region of increased compression is generated, marked with (b) in Figs. 8 and 9. This gas then decompresses and free-falls onto the CPD. However, the majority of the incoming gas is eventually accreted by the planet. This happens either after some residing time in the CPD or by directly flowing over the CPD and colliding with flow from the opposite side of the Hill radius. We note that the global meridional circulation noted in the isothermal runs of Morbidelli et al. (2014) is an azimuthally averaged feature and is therefore not visible in this slice.

loop is in fact a vertical extension of the tidal arms, which must be driven by tidal resonances.

The most prominent vertical structure directly above the planet is the vertical flattening in density at $(x, z) = (0, z)$ which coincides with vertical free-falling gas and the vertical accretion shock seen at $(x, z) = (0, 0.2)$ in the temperature (feature e) in Figs. 8 and 9). This feature is accompanied by a vertical column of high-temperature gas extending up to $(x, z) = (0, 0.5)$ (feature (f) in Figs. 8 and 9). The high-temperature column is a consequence of colliding flows from opposite sides of the planet. This flow collision cancels lateral velocity components and leaves the vertical velocity component non-zero, upon which the gas free-falls vertically onto the planet. This accretion flow is qualitatively similar to what has been postulated and numerically observed in the framework of Bondi accretion (Edgar 2004).

In order to obtain a more quantitative estimate of the effects of those features on local variables, we show the most important variables in vertical 1D profiles in Fig. 11 directly above the planet. From Fig. 11, it is evident that the colliding flows create a negative density gradient $\partial\rho/\partial z$ at $\tilde{z} > 0.35$, but the gas free-falls towards the planet below this level, reversing this gradient. This feature is similar to the 1D accretion shock during runaway gas-accretion, as analysed in Marleau et al. (2017, 2019). We now describe the gas flow that interacts directly with the CPD.

### 4.5. Vertical mass delivery to the CPD I: the free-fall region as an effect of vertically tilted spiral arms

In order to understand the details of the free-fall occurring in the midplane we have to develop an understanding of the behaviour of the mass column directly above it. We see in Sect. 4.2 how the free-fall region is opened as a consequence of overpressure and supersonic erosion in the midplane. In this section, we discuss how the vertical tilt of the spiral arms leads to a cut-off of mass supply to the entire volume of the free-fall region.

Two- or three-dimensional flows that pass-through shocks tilted at an angle with respect to their flow direction behave like oblique shocks, that is, the flow component perpendicular to the shock surface is conserved while that parallel to it experiences

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

**Fig. 11.** Profiles of temperature, vertical velocity, and density as a function of height over the midplane. The letter labels correspond to those used in Fig. 3 and Figs. 8–10. At a height of up to $\tilde{z} = 0.5$ above the planet, we find colliding streams (see also Fig. 10) that originate from overshooting the CPD after originally entering the Hill sphere through the midplane. The resulting compression leads to a slight density bump (f), which is then accreted onto the planet and causes the vertical accretion shock (e).

a discontinuity. As mentioned above, this is often expressed as a vorticity jump, and this vorticity jump depends on the shock obliquity. The vorticity jumps across the spiral arm shock in the midplane of protoplanetary discs (Li et al. 2005, and references therein) is particularly relevant to the physics of gap formation in protoplanetary discs. This is because one can understand the process of gap formation as the effect of the vorticity jump provided by the spiral arm flow patterns.

We find an interesting analogy of this gap opening in the vertical direction. The gas flowing from the CSD is initially stably stratified in the vertical direction. When encountering the tilted shock surface, this gas experiences a vorticity jump that redirects it. The direction into which the gas is redirected depends
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Vertical side view of 3D-integrated streamlines (**left**) going through the circled-star-position in Fig. 3. We note that here we show the $y-z$ plane. Streamlines are coloured according to starting height for purposes of distinguishing them. The circled star (bottom of the plot) is at the same 3D coordinates as before. A flow separation occurs at around $z = 0.41r_H$ dividing the streamlines into strongly shocked upwards and strongly shocked downwards flows. Scanning through horizontally adjacent streamline families reveals that this phenomenon is responsible for leaving the entire column volume behind the shock devoid of streamlines with the ability to replenish the missing mass. Evolution of the free-fall region at the vertical black line (**right**). Density profiles are shown evolving in time at the position indicated by the line in the left figure. The formation timescale of the FFR can be read-off as $\approx 2\Omega^{-1}$. See Fig. 6 for a density map during the formation process.</figcaption></figure>

on the attack angle of the flow towards the shock. This becomes evident from inspecting the results from the 3D integration of a selected family of streamline trajectories in Fig. 12. Those streamlines pass through the spiral arm shock at the position that was marked with a circled star in earlier plots. There, the streamlines that encounter an upwards-tilted shock surface are also tilted upwards, and vice versa. Post shock, the downwards-flowing streamlines are forced to converge through redirection towards the midplane. This causes a region of compression along the flow of streamlines, but most of the streamlines are now directed away from the volume that is geometrically directly behind the spiral arm. We mark this volume simplistically with a black line in Fig. 12.

Once the gas reaches the midplane, a region of high compression is created, which is the ‘foot’ of the shock. At this point, the gas is still ~0.9 $r_H$ away from the planet, but can now fall onto the CPD through the free-fall region.

This vertical flow pattern is repeated for all streamline families that are horizontally adjacent. Hence, the entire columnar volume that forms the free-fall region cannot be replenished with gas by CSD gas due to the vertical tilt of the spiral arm shock (cf. Fig. 12 right). This explains how the lower density state of the free-fall region is maintained, but not what the preconditions are for the vertical tilt of the spiral arm. We investigate this below.

## 4.6. Vertical mass delivery to the CPD II: vertically tilted spiral arms as an effect of competition between ram and static pressure, initiating FFR formation

In order to learn more about the physical conditions behind the spiral arm shock and shed light on the occurrence of the free-fall region, we decided to investigate the pre-shock conditions.

For this purpose, we take a more detailed look at the same family of streamlines as in the last section. We zoom in radially onto the shock, and plot where the transition from locally super-to subsonic gas flow occurs. Because the shock evolves only very

slowly for the first two orbits and then reaches a steady state, the co-rotating frame of the planet and the shock-rest frame coincide. Hence, we analyse the shock data in the co-rotating frame of the planet. The result can be seen in Fig. 13 (left). In this figure, the shock surface can be traced by eye by following the sharp increase and decrease in the vertical velocity component, which are at the points of sharp upwards and downwards turns of the gas, respectively. We refer to this upward and downward transition as the vertical flow separation at a height $z_{sep} = 0.41$. A curious pattern emerges as a function of the vertical direction on the streamlines. The post-shock flow changes its Mach-number from being subsonic to remaining supersonic, at about the height of $z_{sep}$.

In order to analyse this behaviour further, we plot the pre-post-shock jumps for some important shock quantities on the same streamlines, along the shock surface in Fig. 13 (right). First, the jumps of density and temperature are of interest. The bracket $[x]$ denotes the jump ratio $[x] = x_{pre}/x_{post}$. Some physically important values for our parameters are marked in the plot.

In an adiabatic shock with our adiabatic coefficient of $\gamma = 1.4$, the jumps for $[\rho]$ and $[T]$ should approach $[\rho] = 1/4.5$ for $M_{pre} = 4$ (with $[\rho] \to 1/6$ for infinitely strong shocks) and $[T] \approx 1/4$ for $M_{pre} = 4$ (with $[T] \to \infty$ for infinitely strong shocks). From the data it is evident that the behaviour of the spiral arm in the upper layers, for $\tilde{z} > z_{sep}$, is consistent with that of adiabatic shocks.

For $\tilde{z} < z_{sep}$ the shock behaviour changes. An important possibility of modifying the shock jump values is the radiation of the shocked gas. Both sides of the shock radiate energy, and as $[T]$ is unbound, important asymmetries in radiative fluxes can arise. Mihalas & Mihalas (1984, ch. 104) show how the Rankine-Hugoniot jumps for a radiating shock will change compared to the adiabatic case: in a radiating shock, the temperature jump is expected to be more moderate compared to the adiabatic case because both sides of the shock radiate back at each other, while the density jump should be increased.
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Streamlines at the circled star from Fig. 3, with zoom-in on the shock and super-subsonic-transitions (**left**) . The same streamlines are shown as in Fig. 12, with a limited range of azimuthal distances for clarity. We note the asymmetry in the vertical, giving the shock an overall concave shape. Ratios of pre- to post-shock variables are shown for the same vertical position on the shock (**right**). The shock jump values are important for explaining both the shape of the shock and the occurrence of the free-fall region. In the upper parts of the shock, the ram-pressure of the incoming flow is insufficient to overwhelm the shock static pressure of the spiral arm, and the flow is forced to become subsonic (seen where $M_{\text{post}} < 1$). Just below the flow separation at $z_{\text{sep}} = 0.41r_{\text{H}}$, as the density increases, the flow ram pressure becomes sufficient to overcome the static shock pressure. This happens below $z = 0.38r_{\text{H}}$. In this region, the pre-shock supersonic flow also remains supersonic post-shock, albeit slowed down. The flow then continues to advect material from the post-shock volume, leading to the formation of the free-fall region. The region below $z = 0.03r_{\text{H}}$ encounters in its first post-shock cell material coming from above, along the shock curvature. Those colliding flows lead to compressional heating at the “foot” of the shock, increasing the density and temperature contrast. This foot region is important as it delivers mass to the midplane, keeping the midplane density jump alive. Along the entire shock, the adiabatic ratio [$u^2/2 + h$] = 1 holds, revealing that radiative effects play only a minor role for the structure of this shock.</figcaption></figure>

The jump data we see in [$\rho$] and [$T$] for $\tilde{z} < z_{\text{sep}}$ are therefore inconsistent with both adiabatic and radiative shocks, at least using the co-rotating planetary frame as shock-rest frame. Hence, although we see discontinuities in the data, we need to search for another possibility in order to explain this data.

Instead, we show that the ratio of pressures $\zeta$ (as defined in Eq. (3)) explains the data well. Intuitively, it is clear that for $\zeta > 1$ the post-shock gas can remain supersonic. Pushing past the shock should only be possible for incoming gas, if the otherwise unsurmountable pressure barrier of the shock is instead felt as a relatively insignificant “bump on the road”. This is an explanation that seems to agree well with our previous analysis from Sect. 4.2 in relation to Fig. 6 and now in the vertical data in Fig. 13.

The advection of post-shock material leads to the question of why the shock does not disappear altogether and instead appears static. The only explanation we could find for this behaviour was the vertical delivery of mass towards the “foot” of the shock, as evidenced by the additional compressional heating for $z < 0.03r_{\text{H}}$ and shown in Fig. 13. The vertical mass fluxes per streamline are about a factor 30 lower than those advected through the shock in the midplane, but this is the only source of mass available to keep the midplane shock in existence.

In order to spot any anomalies, we plot the total energy jump [$u^2/2 + e + P/\rho$] = [$u^2/2 + h$], where $h$ is the enthalpy, and this yields a constant value of 1 across the whole shock. This confirms that, overall, the shock obeys correct thermodynamics, and the structures we see do not originate in numerical effects.

### 4.7. Vertical mass delivery to the CPD III: spiral arm tilt as an effect of the ram pressure

The dynamic overpressure in the lower parts of the shock, for $\tilde{z} < z_{\text{sep}}$, helps to explain the tilted vertical structure of the spiral

arm structure. The entire region of overpressure is being pushed towards the planet, which distorts the shape of the spiral arm.

We also note that in our other simulation runs, we find agreement with this principle. A convex shape of the spiral arms seems to be a predictor of the existence of a free-fall region, and we now know that the convex shape is tied to the thermodynamics and the intrinsic cooling capability of the envelope. We indicate the qualitative features of the other shocks in Table 2 for simulations C1–C2. Simulations with a higher opacity or a deeper potential generally do not display a FFR due to their increase in pressure support of that region.

### 4.8. Planetary accretion through the CPD, the 3D circulation inside the CPD, and its inner truncation radius

After discussing the mass delivery, we now focus on describing the state of the gas that is entering the CPD from the inner CPD radius before it starts to orbit in the CPD. As noted in Sect. 4.3 and Fig. 7, gas enters mainly through the midplane and is bumped to high values of $v_{\theta}/v_{\text{K}}$ through the action of the CPD spiral arms. This mechanism has more consequences in the vertical, which we explore in this section. Above, we show the general flow structure as a 2D cut in Fig. 10, but for a full understanding of the flows it is necessary to investigate the average circulation inside the Hill sphere, which we show in Fig. 14.

Once passed through the CSD spiral arms, the mass flux enters the Hill sphere predominantly through the midplane, as seen in Fig. 12 (left). Mass fluxes from higher altitudes decrease rapidly by orders of magnitude in strength because the density decreases in altitude and the midplane temperature is low. This results in an integrated per-streamline-flux of 0.1% in the vertical direction, compared to the horizontal.
<table><caption>**Table 2.** Phenomenology of spiral arms, occurrence of free-fall regions, free-fall region evacuation, repeated analysis along the same initial streamline family as for run C1.</caption><thead><tr><th>Run + label</th><th>Spiral arm shape</th><th>Free-fall onto CPD?</th><th>FFR evacuated?</th><th>T<sub>pre</sub></th><th>T<sub>post</sub></th><th>P<sub>tot,pre</sub>/P<sub>stat,post</sub></th></tr></thead><tbody><tr><td>C1 “nominal”</td><td>Concave</td><td>Yes</td><td>Yes</td><td>30</td><td>90</td><td>&gt;2</td></tr><tr><td>C2 “belllin”</td><td>Concave</td><td>Yes</td><td>Yes</td><td>33</td><td>120</td><td>&gt;2</td></tr><tr><td>H1 “high opacity”</td><td>Convex</td><td>No</td><td>No</td><td>70</td><td>90</td><td>1</td></tr><tr><td>H2 “deep”</td><td>Straight</td><td>No</td><td>No</td><td>105</td><td>170</td><td>1</td></tr></tbody></table>

**Notes.** Sample values pre-to-post-shock in order to showcase significant differences. Values are taken above the midplane, to avoid effects of a “foot”-type post-shock region. The pre-shock temperature for the “cold” simulations, C1 and C2 are essentially the disc temperatures. For the hotter simulations H2 and H1, the pre-shock temperatures rise significantly due to radiative precursors.

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Cylindrically averaged vertical structure of flows, density, and temperature in K shown on the contours against planetocentric radius and height. Deviations from a spherically symmetric temperature profile are due to compressional heating for the free-falling and CPD-overshooting midplane flows, and due to efficient radiative cooling in the vertical. Density colours show the toroidal structure of the CPD. Streamline colours indicate the fraction of Keplerian rotation on each flow line, helping to estimate gas residence time on each particular orbit. Streamlines of high Keplerian rotation will therefore orbit more perpendicular to the plot rather than following the streamlines, and vice versa.</figcaption></figure>

The free-falling gas coming from the CSD and the spiral arms enters the CPD after being slowed down in the midplane accretion shock, which shows a concave shape, the opposite of the spiral arms. This can be seen in Fig. 14, together with the associated density and temperature structures of the CPD. There, vertical kinks in the temperature contours indicate changes in the compressional heating of flows, while the density structure of the CPD is relatively simple, at a high aspect ratio of $H/r \approx 0.1-0.2$. Streamlines that pass the CPD boundary at around $r \approx 0.4-0.5r_H$ exit again in the vertical direction, but further inward, starting at $(\tilde{r}, \tilde{z}) \approx (0.25, 0)$. This exit is a comparatively slow process, as the gas is redirected on circular orbits with high $v_\theta/v_k$ after the CPD accretion shock and only slowly spirals upwards and outwards from the planet for ~10–20 CPD orbits. This process terminates at the height where the CPD spirals end (previously feature (d) in Fig. 9 and here at $(\tilde{r}, \tilde{z}) \approx (0.4, 0.3)$). This is indicative that this slow upward and outward migration of streamlines seems to be related to the kicks by the CPD spirals at every orbit, similar to the boosting that has been seen in the midplane (compare this to Fig. 7) but with a vertical component added to it. After rising to the height of CSD spiral termination, the gas moves towards the planet much faster.

This gives the gas enough time to cool upon reaching higher into the CSD atmosphere. Viscosity allows the streamlines drift inwards slowly until the streamlines reach the evacuated region directly above the planet and then fall towards it. Of those

streamlines that reach the inner ~0.2$r_H$ a small number of them goes into a circulatory motion caught between $r_s = 0.1r_H$ and $0.3r_H$, corresponding to the circulation seen earlier in the midplane in Fig. 7.

Mass accretion through the midplane therefore does not seem to deliver mass very efficiently to the CPD: the gas is merely passing through, albeit in a fairly complex manner, until it finally lands inside the smoothing length of the planet. Except for the small region between $r_s = 0.1r_H$ and $0.3r_H$, which can compensate for lower Keplerian rotation with radial–vertical circulation, we do not find any region in the CPD which could possibly permit long residence times for dust particles. This indicates that our simulation results pose a too early stage of the CPD life in order to form moons.

## 5. Discussion and support simulations: comparison with entire data set, CPD masses, and planetary accretion rates

The previous section describes the physical properties of the gas flow, density, and temperature for the nominal simulation run. We now begin to introduce and compare the simulation runs for a number of important parameter variations. On one hand, for the discussion of CPD properties, we limit ourselves to comparing runs C2-H2 to C1, while on the other hand, in order to place our
results into a larger context of growing protoplanets, we discuss mass accretion rates into all planetary envelopes and include runs m1-m5 in the discussion.

## 5.1. Comparison to other simulation runs and robustness of C1 results

In order to simplify the comparison of run C1 to our other Jupiter-mass simulation runs, we avoid plotting the full complexity of data shown previously. We only plot 1D profiles in the midplane and the vertical as this is sufficient to comment on the structure of the envelopes of the simulations C2–H2 in Fig. 15.

The constant opacity of run H1 is hundred times higher than that of C1. This parameter change is motivated by the possibility of dust-enriched, ISM-like material flowing into the feeding region of the giant planet being distributed into the gas and hence increasing the ambient dust opacity. It also serves as an important comparison to earlier work using such high opacities (Szulágyi & Mordasini 2017; Lambrechts et al. 2019). The physical extension of the envelopes in H1 is large enough to preclude large free-fall distances, and additionally the temperature is everywhere high enough so that none of the infalling gas becomes supersonic. Hence, the paucity of physical features in H1 (see the top row of Fig. 15) is consistent with those earlier works, while the Hill sphere is still accreting vigorously at a rate of $\dot{m}_{H1} = 1.0 \times 10^{-2} m_{\oplus} \text{yr}^{-1}$. The conclusion that those features are physical rather than numerical is supported by the fact that we ran H1 with the same numerical parameters as C1, the only difference being the opacity.

The potential of run H2 is four times deeper than that of C1 (and H2 has correspondingly higher numerical resolution). This potential depth of $\tilde{r}_s = 0.025$ corresponds to a physical size of $r_s = 12 R_{\text{Jup}}$. This physical size, while still being too large for a young protoplanet, is an improvement over our nominally used potential depths and those of our previous works, but limits the simulation speed significantly. The deeper potential causes strong initial compression of the infalling gas and hence high temperatures throughout the whole envelope. The central temperatures reach about 13 000 K, which is very high, but this may be an artifact of our constant $\gamma$ treatment that ignores the dissociation of the hydrogen molecule. The larger extent of the planetary boundary condition inwards serves as a check for the robustness of our inner CPD boundary.

As can be seen in Fig. 15 (middle row) the emergence of a free-fall region is suppressed when deepening the potential in H2. We find that this is consistent with our previous analysis of spiral shock shape and pressure ratio $\zeta$. The intense radiation field released by the deep potential is shining on the spiral arm and stabilises it by increasing the static pressure. Data on the spiral arm phenomenology for runs C1 and C2–H2 are summarised in Table 2.

The rotational state of the H2 envelope is only weakly impacted by the elevated temperature. As can be seen in Fig. 16 (left), the envelope achieves above 80% Keplerian rotation. Hence, contrary to Szulágyi et al. (2016) we obtain Keplerian rotation values that are as high as those in the cold simulation of these latter authors, but with the planetary temperatures of their hot simulation. This is because the lower opacities compared to Szulágyi et al. (2016) allow for higher temperature gradients between the deep planetary envelope and $r_H$, meaning that the planet can be very hot but still allow for the CPD to cool down. It seems that the central planetary temperature does not determine

whether a CPD or envelope forms, but rather how far inwards it extends.

The CPD also possesses a flow separation between inner and outer envelope at $r \approx 0.1 r_s$, at the same radius as run C1. This demonstrates the robustness of this inner CPD edge against our numerical parameters. However, the global flow inside the Hill sphere has changed significantly, as is evident from Fig. 16 (right).

This global flow pattern in H2 follows a global circulation pattern, bearing resemblance to that seen in Tanigawa et al. (2012) (obtained with locally isothermal simulations), and is reasonably similar to that seen in H1, which is why we do not show the latter simulation separately. The reason for the similarity does not seem to be the vertically quasi-isothermal nature of the H2 CPD, as the vertical temperature profile is similar in C1, but rather the radial structure of the pressure gradients. Interestingly, the accretion rates do not differ strongly when comparing C1 and H2. The accretion rate onto the H2 envelope is $\dot{m}_{H2} = 1.8 \times 10^{-2} m_{\oplus} \text{yr}^{-1}$ whereas the accretion rate onto the C1 envelope is $\dot{m}_{C1} = 3.5 \times 10^{-2} m_{\oplus} \text{yr}^{-1}$.

Finally, run C2 uses the more realistic opacities of Bell & Lin (1994) instead of constant opacities. This set of opacities features two important opacity transitions, which correspond to the sublimation temperatures and densities of water and silicates. This run has a long burn-in time of the simulation because of the non-linear time evolution of the envelope opacities, but eventually finds a quasi-steady state. While the exact values of density and temperature in C2 differ slightly compared to C1, we recover all simulation features of C1 also with the non-constant Bell & Lin opacities. This is an important data point showing the validity of our previous analysis.

We feel it necessary to stress that the global flow patterns inside the Hill spheres of our Jupiter-mass planets seem to divide along two different kinds: C1 and C2 both feature a free-fall region due to their low envelope temperatures and their flows are very similar and complex. Simulations H2 and H1, being “hotter”, feature a simple global circulation accreting significant mass through the vertical direction, albeit with the same inner CPD cut-off as the ‘cold’ circulations.

For both types of circulation, the accretion rates remain high and this is discussed further in the following section. Most of the mass (≈99%) entering the Hill sphere enters through the midplane and is then either processed through the CPDs, after residence times of several tens of CPD orbits, or evades the CPD and falls directly onto the planet at a vertical angle of ≈45°.

Those simulation results demonstrate that a low mean opacity is key for the development of circumplanetary discs. In star-forming environments this could be achieved by having a reduced total dust mass of micron-sized grains, which dominate cooling rates if present. Another possibility to obtain lowered mean opacities is via growing the same dust mass to millimetre(mm) sizes on average. Evidence for dust growth in prestellar cores (Chacón-Tanarro et al. 2019) has been reported, and even the presence of mm-sized grains has been inferred indirectly in protoplanetary discs (Harsono et al. 2018). Therefore, our choice of mean opacity values is well supported by observational evidence.

## 5.2. Gas accretion into CPD, CPD masses, and radial mass profiles

While the planetary accretion rates are vigorous for the Jupiter-mass planets, it is interesting to note that our CPDs also accrete,
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Overview of slices through flow and structure variables, identical to Figs. 4 and 11, only now for the runs H1 ("high opacity", top), H2 ("deep", middle), and C2 ("Bell & Lin", bottom) at orbit 5 in the midplane (left column) and the vertical values (right column). The letter labels indicate the same physical features as in previous figures: CSD spiral arm shock (a), same volume as the free-fall region in run C1 for comparison (b), shock due to free-fall of gas from the spiral arms towards the CPD (c), temperature bump due to the CPD spiral arms (d), vertical accretion shock (e), and accretion column from colliding flows (f). It becomes evident that the envelope with high opacity suppresses structural and flow features, consistent with the literature (e.g. Lambrechts et al. 2019). A deeper potential raises temperature even above those of the low-opacity run, but only suppresses the formation of a free-fall region. The Bell & Lin opacities recover all features we found previously in the nominal run.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Midplane rotational structure for the run H2 (**left**), with $|v|/v_{Kep}$ as contours. The potential of run H2 is four times deeper than that of run C1, which is shown in Fig. 7. The gravitational smoothing length is marked with the blue dashed circle. It becomes evident that the deeper potential does not extend the CPD necessarily much closer to the planet, rather the contrary. As pressure support increases, the maximum rotational support drops below 90%. The physical extent of the strongly rotating gas remains similar to C1. Cylindrically averaged structures are shown with the temperature as contours (**right**). The observed flows are similar to the isothermal results presented in Tanigawa et al. (2012) and consistent with previous high-opacity results from Lambrechts et al. (2019).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Circumplanetary disc masses for the Jupiter-mass planets (**left**) seem to remain stable after the initial two orbits of ramp-up time of $r_s$. The CPDs, with the exception of C1, start accreting with gas accretion rates in the range of $1-4 \times 10^{-5} m_{\oplus} \text{yr}^{-1}$, corresponding to growth times for the Galilean satellites between 80 kyr and 0.32 Myr. Simulations are not run for the same number of orbits due to runtime restrictions. Mass accretion rates onto the planets as a function of planet mass and opacity (**right**). We note similarities to Lambrechts et al. (2019), which used Bell & Lin opacities for a similar mass range. Earlier results from Ayliffe & Bate (2009b) using SPH simulations and interstellar grain opacities (IGO) agree with ours in the mid-mass range, but our low masses accrete visibly more and the high masses do not seem to converge on a single disc-limited accretion rate. One could speculate that this is related to the strong disc inflow that triggers the CSD spiral arm tilt, effectively increasing the mass inflow.</figcaption></figure>

which can be seen in Fig. 17 where we plot the CPD masses, defined as the mass between the shells of $0.1 = \tilde{r}_s < \tilde{r} < 0.5$. One conclusion from the sections above was that the CPD in runs C1 and C2 is mainly a complicated pass-through with very low permanent storage; it now becomes evident that some mass is accumulated in this region.

The only simulation that does not seem to accrete mass onto the CPD is C1. This is interesting, as the flow pattern seen in this run is very similar to that of C2. The significant mass noise in this region makes it difficult to decipher whether C1 is slowly accreting or not at all.

The radial distributions of CPD surface densities, parameterised as $\Sigma(r) = \Sigma_0 r^a$, seem to be impacted only weakly by opacity effects, but are dominated by the planetary potential depth. We use $\Sigma_0 = \Sigma(0.2 r_H) \approx 0.7 \text{ g cm}^{-2}$, a value which is approximately the same for all runs at this radius. We find relatively steep power-law functions for the low-opacity runs C1 and C2 of $a = -2.5 \pm 0.2$ in the range $0.4 < \tilde{r} < 1.0$, which reflect the influence of the free-fall region on the mass profile, and a significantly flatter $a = -1.0 \pm 0.2$ in the range $0.1 < \tilde{r} < 0.3$, with some more fluctuations for C2 due to opacity changes; Fujii et al. (2017) found similar deviations from the power-law
behaviour when using non-constant opacities. A slightly shallower distribution is found for run H1 with $a = -2.0 \pm 0.2$ over the whole envelope. However, the latter fits are only extended over the relatively modest CPD region in those runs. A more realistic power-law is probably given by the CPD for run H2, for which we find a significantly flatter power-law with $a = -1.0 \pm 0.2$ for $0.1 < \tilde{r} < 1.0$ and in the pressure-influenced region $a = -2.0 \pm 0.2$ for $0.025 < \tilde{r} < 0.1$.

The gas residence times in all CPDs are fairly low, and therefore the mass profiles are a direct consequence of the mode of radial mass transport that each individual CPD finds on this short timescale. A future set of simulations will be dedicated to studying deep potentials like that in H2, for a wider range of parameters. However, it is encouraging that the value for *a* that we find in the CPD region of the nominal simulation run is extended over a much larger range in radii in run H2, suggesting that the true value for *a* in our simulation settings is $a \approx -1$.

The runs with a surface density power law of $a \approx -2$ have most of their mass stored in the inner CPD, which is unresolved in those simulations. Run H2, which has a more realistic potential depth but is very hot, with $a \approx -1$, has a CPD mass of $m \sim \int_{r_1}^{r_2} dr r \Sigma(r) \propto \Sigma_0(r_2 - r_1)$. Therefore, the maximum CPD mass for this mass profile is $\Sigma_0 r_2 = \Sigma_0 r_H$. The value of $\Sigma_0$ is set by the mass flowing through the CPD, but it is conceivable that once the CPD cools down, the lifetime of the gas increases significantly in the CPD which would stabilise the mass profile and allow for higher CPD accretion rates.

Our measured CPD masses are very low compared to what is commonly assumed in models for the formation of moon systems (Canup & Ward 2006; Cilibrasi et al. 2018; Ronnet & Johansen 2020). The ratio between our CPD masses and the central planetary mass is of the order of $10^{-5}$. Isella et al. (2019) reported a CPD mass for PDS 70c that is ≈100 times larger than our values. It is possible that this is a consequence of the surface density profile and the inner CPD truncation radius. However, with the measured values for the CPD accretion rates (see Fig. 17) the needed dust mass can be delivered in ≈10 kyr, while for the Galilean satellites between 80 kyr and 0.32 Myr would be required.

### 5.3. Gas accretion rates into the planetary envelopes

After discussing the CPD accretion rates, we now analyse the envelope accretion rates. We list the measured CPD masses for the Jupiter-mass planets and the average accretion rates over the simulation time into the Hill sphere for all our planets in Fig. 17. Our results show that for well-resolved simulations there is significant gas accretion for all the considered opacities and for all planetary masses. This importantly implies that there is no problem in 3D gas accretion simulations once they are well resolved, unlike previously stated in Szulágyi et al. (2014).

The accretion rate that we find for the $240 m_⊕$ planets for both opacities is slightly higher than for the $1 m_J$ planet, as can be seen in Fig. 17. From this we conclude that the Jupiter-mass planets are in a state of disc-limited accretion for the low-opacity case, that is, the accretion rate is no longer limited by the Kelvin-Helmholtz contraction but rather by what the protoplanetary disc can provide through the Hill sphere. We refer to Lambrechts et al. (2019) for further analysis of disc-limited versus cooling-limited accretion. This is classically understood to be a process in which the Kelvin-Helmholtz contraction of the envelope material is faster than any replenishment from the disc. This results in the planet acting to accrete all the gas supplied through the

Hill sphere. The data points presented here differ quantitatively from the classic results in Ayliffe & Bate (2009b) mainly in that our accretion curve is flatter towards the low masses and we do not see a strong flattening of accretion rates for the high masses. Their results were obtained using an SPH code, and it is possible that the differences from our *ṁ(m)*-curves result from differences in implementing radiation transport in Eulerian versus Lagrangian codes.

In the context of the work of Cimerman et al. (2017), it is perhaps surprising that we do not see any effects of entropy recycling at the low-mass end in our results, which would stop the gas accretion process. The simulations of these latter authors, while well-resolved, operate at lower planet masses (their maximum planet mass is $5 m_⊕$) and at much higher optical depths as well as disc masses (their $\rho_{\text{mid}} \approx 6 \times 10^{-6}$ g cm$^{-3}$, while we have $\rho_{\text{mid}} \approx 4 \times 10^{-11}$ g cm$^{-3}$ for the unperturbed discs of the $20 m_⊕$ planets). It is because of this main difference that we do not see any important recycling effects for our low-mass planets.

## 6. Summary

Here, we present a set of 3D radiation hydrodynamics simulations that were created and analysed in order to learn about the occurrence of circumplanetary discs and mass delivery onto them.

1. Initially, we asked at which planet masses CPDs would occur and whether or not opacity is an important factor. We find that there is a significant shift in the rotation profiles of envelopes between 120 $m_⊕$ and 240 $m_⊕$, indicating the takeover of rotational support from pressure support. This indicates that CPDs form only late in the envelope accretion process of a gas giant. The residence time of gas in the CPDs is of the order of a few tens of CPD orbits; hence the rotation of the CPD gas is a direct result of newly accreted gas and not inherited from the initial conditions. We find that low opacity significantly helps flatten envelopes and reduce pressure support. Surface density profiles in regions with significant rotational support stabilise to $\Sigma \sim r^{-1}$.

2. Jupiter-mass planets accrete at disc-limited accretion rates of $10^{-2} m_⊕$ yr$^{-1}$. Mass enters the planetary Hill sphere predominantly through the midplane, and after being processed in the CPD ends up on the planet. This remains true even after extending the potential depth of our planets by a factor of four. Circumplanetary discs seem to accrete with values of around $1-4 \times 10^{-5} m_⊕$ yr$^{-1}$, but those values are difficult to determine and there is a possibility that calculations are susceptible to numerical effects.

3. The flows inside the CPD are influenced by the occurrence of a separate set of CPD spiral arms, as opposed to the classical CSD spiral arms which are previously identified in the literature as tidal arms (Zhu et al. 2016). Those CPD spiral arms cause gas to slowly rise above the midplane while it orbits, and finally to spiral into the planet at high altitude. The CPD spiral arms therefore play a role in shaping the CPD structure and channelling the accreted gas into the planet.

4. In our investigation of a low-opacity Jupiter (run C1), we found a complete planetary detachment from the parent CSD expressed as a region of pressure gradient inversion. We find a new circulation regime inside the Hill sphere. In this case, the spiral arms feature sufficiently low static pressure so as to be overcome by the ram pressure of the accretion flow. This ram pressure evacuates a region between the spiral arm and the CSD. This newly evacuated region facilitates free-fall
onto the CPD, which causes a weak accretion shock in the midplane, additional to a pre-existing vertical accretion shock. The spiral arms are tilted through the overpressure, which can be used as a phenomenological identification mark for this effect. We note that this effect depends on the efficient radiative cooling of the spiral arm and functions at low envelope temperatures.

5. Using low Bell & Lin (1994) opacities (run C2) instead of constant ones reproduces a very similar simulation outcome to the one seen in C1.

6. High opacity (run H2) reproduces other work from the literature, with this envelope being an almost featureless blob, with a large-scale circulation inside the Hill sphere. Nonetheless the accretion rate for this run is on the order of $10^{-2} m_{\oplus} \text{ yr}^{-1}$ into the smoothing length region.

7. A simulation run with four-times-deeper potential than that of C1, but otherwise identical parameters, was investigated (run H1). We find very hot temperatures ($T_c \approx 13\,000\,\text{K}$), but significant rotation in the outer CPD ($v_\theta/v_K \approx 0.8$) and an identical extent of the CPD to inner radii as in C1. This indicates the robustness of our findings for the CPD structure. Furthermore, a lowered opacity, which can be obtained under a number of realistic circumstellar conditions (Ossenkopf & Henning 1994) from various dust populations, seems to be paramount for forming CPDs.

A major limitation of our model is nevertheless the limited numerical resolution of the inner regions of the Hill sphere. We observe that increasing the resolution decreases the luminosity of the planet and significantly affects the temperature and structure of the CPD. Therefore, future investigations of CPDs should be performed at higher resolutions in order to identify resolution convergence criteria regarding their structure.

The constant turbulent viscosity of $\nu \approx 10^{15} \text{ cm}^2 \text{ s}^{-1}$ used here, corresponding to $\alpha \approx 10^{-2}$, is relatively high in light of recent observational evidence (Pinte et al. 2016; Dullemond et al. 2018). A re-examination of our results for lowered viscosities is therefore being carried out. However, due to the significant run times of the gap formation at low viscosities in 3D RHD, we were not able to include the results here.

*Acknowledgements.* We thank the anonymous referee for improving the quality of our manuscript and discussions of shock physics. M.S. was supported by a project grant from the Swedish Research Council (grant number 2014-5775). A.J. wants to thank the support by the Knut and Alice Wallenberg Foundation (grant number 2017.0287), the European Research Council (ERC Consolidator Grant 724 687-PLANETESYS) and VR grant 2018-04867. M.L. thanks the Knut and Alice Wallenberg Foundation under grant 2017.0287. B.B. thanks the European Research Council (ERC Starting Grant 757 448-PAMDORA) for their financial support. All the simulations presented in this work were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lundarc in Lund University, Sweden, and the entire team is grateful for being supported with their expertise.

## References

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

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

Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987

Béthune, W. 2019, MNRAS, 490, 3144

Birnstiel, T., Andrews, S. M., & Ercolano, B. 2012, A&A, 544, A79

Bitsch, B., Morbidelli, A., Lega, E., & Crida, A. 2014, A&A, 564, A135

Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859

Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834

Chacón-Tanarro, A., Pineda, J. E., Caselli, P., et al. 2019, A&A, 623, A118

Cilibrasi, M., Szulágyi, J., Mayer, L., et al. 2018, MNRAS, 480, 4355

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

Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2

D'Angelo, G., & Lubow, S. H. 2008, ApJ, 685, 560

D'Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540

Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46

Edgar, R. 2004, New Astron. Rev., 48, 843

Fujii, Y. I., Okuzumi, S., & Inutsuka, S.-i. 2011, ApJ, 743, 53

Fujii, Y. I., Okuzumi, S., Tanigawa, T., & Inutsuka, S.-i. 2014, ApJ, 785, 101

Fujii, Y. I., Kobayashi, H., Takahashi, S. Z., & Gressel, O. 2017, AJ, 153, 194

Goldreich, P., & Tremaine, S. 1978, ApJ, 222, 850

Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425

Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59

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

Harsono, D., Bjerkeli, P., van der Wiel, M. H. D., et al. 2018, Nat. Astron., 2, 646

Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25

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

Kevlahan, N. K.-R. 1997, J. Fluid Mech., 341, 371

Klahr, H., & Kley, W. 2006, A&A, 445, 747

Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971

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

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

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

Lega, E., Crida, A., Bitsch, B., & Morbidelli, A. 2014, MNRAS, 440, 683

Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

Li, H., Li, S., Koller, J., et al. 2005, ApJ, 624, 1003

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

Marleau, G.-D., Mordasini, C., & Kuiper, R. 2019, ApJ, 881, 144

Masset, F. 2000, A&AS, 141, 165

Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (New York: Dover Publications)

Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266

Mordasini, C. 2014, A&A, 572, A118

Ormel, C. W., Kuiper, R., & Shi, J.-M. 2015, MNRAS, 446, 1026

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Owen, J. E., & Menou, K. 2016, ApJ, 819, L14

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Ronnet, T., & Johansen, A. 2020, A&A, 633, A93

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

Shibaike, Y., Ormel, C. W., Ida, S., Okuzumi, S., & Sasaki, T. 2019, ApJ, 885, 79

Szulágyi, J., & Mordasini, C. 2017, MNRAS, 465, L64

Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65

Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853

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

Zhu, Z., Ju, W., & Stone, J. M. 2016, ApJ, 832, 193
# Appendix A: The origin of the spiral arm tilt

The vertically tilted spiral arms documented here are the centre point of the complex physics which we analyse in this work. As their structure seems to be related to the cooling properties of the envelope and the disc, it is natural to pursue an explanation of the spiral arm tilt in the vertical temperature structure.

Therefore, in Fig. A.1 we document the unperturbed temperature profiles in the disc at the position of the future planet and at the positions of the future gap edges. This simulation is run in 2D for 400 orbits in order to quickly find a radial–vertical equilibrium. The temperature profiles are Gaussian, as expected for a quasi-isothermal disc.

To put those into context, we also plot the temperatures for our nominal high-resolution simulation run C1 (dashed lines in Fig. A.1). In view of the vertical temperature inversion, we initially assumed to have found the reason for the static pressure argument from Fig. 13. However, a control simulation with larger radial box size (dotted lines in Fig. A.1) evolved to a different temperature structure with a planet due to a change in the flow in the upper disc layers (compressional heat plays an important role at those low temperatures). This temperature structure shows no temperature inversion; however, the spiral arm tilt, with associated free-fall region, was still evident in the control simulation.

This led us to reject the hypothesis that the spiral arm tilt is caused directly by the unperturbed disc vertical temperature gradients, and we presume that it must be the result of the self-consistently computed temperature solution in the envelope. We further note for possible future studies, that the spiral arm tilt is already evident in the gap-forming run. This should make it possible to quickly find good candidate envelopes if one wants to study the spiral arm tilt without having the need to scan parameter space in expensive high-resolution simulations.

<figure><img src="image_19.png" /><figcaption>**Fig. A.1.** Vertical, unperturbed temperature profile after 400 orbits before the planet is injected (solid red, blue, and black lines) plotted against the distance to the midplane in units of the future planet. Radial positions of the profiles are given in the plot, the azimuthal position is always $\theta=0$. The abrupt temperature decrease at the upper simulation boundary is due to the 3K boundary condition, but does not influence the simulation result because the disc there is optically thin. Dashed and dotted lines are after five orbits of the simulation C1 and a control simulation with a larger radial box size, respectively. Blue for the $r=0$ in the high-resolution runs is not plotted, as it appears in Figs. 11 and 15.</figcaption></figure>