© ESO 2020

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 bottlenecks 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 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 circumplanetary 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⁻¹. Viscous momentum transport is computed through the full hydrodynamic stress tensor multiplied with a constant physical viscosity of ν = 10⁻¹⁵ cm² s⁻¹ that corresponds roughly to α ≈ 10⁻² 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 ≈ 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 (κ = 1.0 cm² g⁻¹) 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_{\text{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}}_{\text{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}}_{\text{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 κ or those from Bell & Lin (1994) scaled by a factor, which we list as ɛ 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² to 10⁴ 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, θ, ϕ]. We also use the planetocentric coordinates [x, y, z] = [(r − 1) cosθcosϕ, r sinθcosϕ, r sinϕ], 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 θ ∈ [−π, π], $\varphi \in {[81,90]}^{°}$ 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 ϕ = 90° for all variables. Reflective boundaries are used for hydrodynamic variables at ϕ = 81°, and open boundaries for the radiative energy at ϕ = 81°. 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 ∈ [0.7, 1.3] × 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 ∈ [0.4, 1.6] × 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Ξ ∝1 − g(Ξ − Ξ₀), where g(Ξ) is a Gaussian function, Ξ represents any of the three spatial directions r, θ, ϕ in our spherical coordinate system, and Ξ₀ 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Ξ_min, that is, the element with the best resolution directly adjacent to the planet and maximum of g(Ξ), 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 × 10⁻⁴ (or 7.5 r_Jup) in size and our simulation run H2 with the deep potential is of dx_min = 1.7 × 10⁻⁴ (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 support, the aspect ratio can only be recast into the form $$\frac{H}{r}={\left(\frac{c_{\text{s}}}{v_{\text{k}}}\right)}^2,$$(1)

which is not v_θ∕v_k. Hence, in order to measure the flatness of a disc we define the flatness parameter at distance d through the cylindrically averaged 2Ddensity $\overline{\rho }(r,z)\text{}:$ $$\xi (d)\equiv \frac{\overline{\rho }(0,d)}{\overline{\rho }(d,0)}.$$(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 ξ additionally to v_θ ∕v_k to assess the properties of planetary envelopes with increasing mass for the simulation sets m1–m5, C1, and H1.

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 mp ≥ 240 m⊕. The strict transition into the Keplerian shear of the CSD happens at vθ∕vKepler=−1, as can be seen from the slope turnover in the rotation curves. The vertical black dashed curve denotes r = 0.3rH, where rotation and density asymmetry values are measured and plotted in Fig. 2.

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 allthe 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, θ_p, ϕ_p) = (1, 0, π∕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_⊕ and 60 m_⊕ 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_⊕ to 120 m_⊕ there is only a relatively weak evolution in the envelopes. After a mass-doubling from 120 m_⊕ towards 240 m_⊕, a transition in the shapeof the rotational profile occurs, developing a “pedestial” consisting of relatively high Keplerian rotation in the CPD region for the 240 m_⊕ and 360 m_⊕ 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.

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_θ ∕v_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 κ = 0.01 cm² g⁻¹ 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).

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θ ∕vK 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⊕ 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.

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.

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 midplane 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 temperatureas 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 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 ≈ 30 K.

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.

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 Ω⁻¹. 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 supersonicafter 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}}},$$(3)

where u_⊥ corresponds to the streamline-velocity component perpendicular to the shock surface. We stress at this point that the nominator and denominator of ζ are generally independent of each other due to the supersonic nature of the inflow from the circumstellar disc. We use the definition of ζ 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 ζ < 1. Those that remain supersonic post-shock coincide with the area at the shock where ζ > 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 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² 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.

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 Ω 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Ω−1 showcase the slow evolution of the spiral arms between those times. The black star denotes the same position as in Fig. 3.

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 Ω−1. 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 themiddle 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.

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 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 long time inthe 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.

Fig. 7 Midplane rotational structure for the run C1 with |v|∕vKep 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|∕vKep. This resultsin 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.

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 arelabelled 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).

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.

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 previously 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 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 ∂ρ∕∂z at $\tilde{z}>0.35$, but the gas free-fallstowards 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.

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 eventuallyaccreted 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.

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 overpressureand 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 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 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.

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 isthen accreted onto the planet and causes the vertical accretion shock (e).

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.41rH 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 ≈2 Ω−1. See Fig. 6 for a density map during the formation process.

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- to 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 γ = 1.4, the jumps for [ρ] and [T] should approach [ρ] = 1∕4.5 for M_pre = 4 (with [ρ] → 1∕6 for infinitely strong shocks) and [T] ≈ 1∕4 for M_pre = 4 (with [T] →∞ 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_{\text{sep}}$, is consistent with that of adiabatic shocks.

For $\tilde{z}<z_{\text{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.

The jump data we see in [ρ] 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 ζ (as defined in Eq. (3)) explains the data well. Intuitively, it is clear that for ζ > 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_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 + e + P∕ρ] = [u²∕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.

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 Mpost < 1). Just below the flow separation at zsep = 0.41rH, as the density increases, the flow ram pressure becomes sufficient to overcome the static shock pressure. This happens below z = 0.38rH. 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.03rH 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 [u2 ∕2 +h] = 1 holds, revealing that radiative effects play only a minor role for the structure of this shock.

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.3and Fig. 7, gas enters mainly through the midplane and is bumped to high values of v_θ ∕v_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.

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 ≈ 0.1–0.2. Streamlines that pass the CPD boundary at around r ≈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_θ ∕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.2r_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.

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 flowline, 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.

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 ṁ_H1 = 1.0 × 10⁻² m_⊕ yr⁻¹. 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}}_{\text{s}}=0.025$ corresponds to a physical size of r_s = 12 R_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 γ 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 ζ. 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 ≈ 0.1r_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 ṁ_H2 = 1.8 × 10⁻² m_⊕  yr⁻¹ whereas the accretion rate onto the C1 envelope is ṁ_C1 = 3.5 × 10⁻² m_⊕  yr⁻¹.

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 theHill 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-formingenvironments 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, 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 Σ(r) = Σ₀r^a, seem to be impacted only weakly by opacity effects, but are dominated by the planetary potential depth. We use Σ₀ = Σ(0.2 r_H) ≈ 0.7 g cm⁻², 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 ± 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 ± 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 ± 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 ± 0.2 for $0.1<\tilde{r}<1.0$ and in the pressure-influenced region a = −2.0 ± 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 ≈−1.

The runs with a surface density power law of a ≈−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 ≈−1, has a CPD mass of $m~{\displaystyle {\int }_{r_1}^{r_2}\text{d}}rr\Sigma (r)\propto {\Sigma }_0(r_2-r_1)$. Therefore, the maximum CPD mass for this mass profile is Σ₀ r₂ = Σ₀ r_H. The value of Σ₀ 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⁻⁵. 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.

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.

Fig. 16 Midplane rotational structure for the run H2 (left), with |v|∕vKep 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).

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 rs. The CPDs, with the exception of C1, start accreting with gas accretion rates in the range of 1–4 × 10−5 m⊕ 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 inthe 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.

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 ρ_mid ≈ 6 × 10⁻⁶ g cm⁻³, while we have ρ_mid ≈ 4 × 10⁻¹¹ g cm⁻³ for the unperturbeddiscs 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 Σ ~ r⁻¹.

• 2.

  Jupiter-mass planets accrete at disc-limited accretion rates of 10⁻² m_⊕ yr⁻¹. 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 × 10⁻⁵ m_⊕ yr⁻¹, 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⁻² m_⊕ yr⁻¹ 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 ≈ 13 000 K), but significant rotation in the outer CPD (v_θ∕v_K ≈ 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 ν ≈ 10¹⁵ cm² s⁻¹ used here, corresponding to α ≈ 10⁻², 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 Lunarc in Lund University, Sweden, and the entire team is grateful for being supported with their expertise.

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 temperatureinversion; however, the spiral arm tilt, with associated free-fall region, was still evident in the control simulation.

This ledus 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.

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 θ = 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.

References

All Tables

All Figures

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 mp ≥ 240 m⊕. The strict transition into the Keplerian shear of the CSD happens at vθ∕vKepler=−1, as can be seen from the slope turnover in the rotation curves. The vertical black dashed curve denotes r = 0.3rH, where rotation and density asymmetry values are measured and plotted in Fig. 2.

In the text

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θ ∕vK 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⊕ 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.

In the text

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.

In the text

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.

In the text

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 Ω 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Ω−1 showcase the slow evolution of the spiral arms between those times. The black star denotes the same position as in Fig. 3.

In the text

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 Ω−1. 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 themiddle 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.

In the text

Fig. 7 Midplane rotational structure for the run C1 with |v|∕vKep 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|∕vKep. This resultsin 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.

In the text

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 arelabelled 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).

In the text

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.

In the text

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 eventuallyaccreted 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.

In the text

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 isthen accreted onto the planet and causes the vertical accretion shock (e).

In the text

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.41rH 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 ≈2 Ω−1. See Fig. 6 for a density map during the formation process.

In the text

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 Mpost < 1). Just below the flow separation at zsep = 0.41rH, as the density increases, the flow ram pressure becomes sufficient to overcome the static shock pressure. This happens below z = 0.38rH. 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.03rH 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 [u2 ∕2 +h] = 1 holds, revealing that radiative effects play only a minor role for the structure of this shock.

In the text

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 flowline, 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.

In the text

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.

In the text

Fig. 16 Midplane rotational structure for the run H2 (left), with |v|∕vKep 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).

In the text

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 rs. The CPDs, with the exception of C1, start accreting with gas accretion rates in the range of 1–4 × 10−5 m⊕ 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 inthe 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.

In the text

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 θ = 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.

In the text