Issue 
A&A
Volume 690, October 2024



Article Number  A183  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202450899  
Published online  08 October 2024 
Gas dynamics around a Jupitermass planet
I. Influence of protoplanetary disk properties
^{1}
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
Laboratoire Lagrange,
France
^{2}
Univ. Grenoble Alpes, CNRS, IPAG,
38000
Grenoble,
France
^{3}
UniversitätsSternwarte München, LudwigMaximiliansUniversität,
Scheinerstr. 1,
81679
München,
Germany
^{4}
Collège de France,
11 Pl. Berthelot,
75005
Paris,
France
^{5}
Imperial College London,
Kensington Lane,
London
W1J 0BQ,
UK
^{6}
Department of Earth, Planetary, and Space Sciences, University of California,
Los Angeles,
CA90095,
USA
^{7}
Center for Star and Planet Formation, Globe Institute, University of Copenhagen,
ØsterVoldgade 57,
1350
Copenhagen,
Denmark
^{★} Corresponding author; email: elena@oca.eu
Received:
28
May
2024
Accepted:
3
August
2024
Context. Giant planets grow and acquire their gas envelope during the disk phase. At the time of the discovery of giant planets in their host disk, it is important to understand the interplay between the host disk and the envelope and circumplanetary disk properties of the planet.
Aims. Our aim is to investigate the dynamical and physical structure of the gas in the vicinity of a Jupitermass planet and study how protoplanetary disk properties, such as disk mass and viscosity, determine the planetary system as well as the accretion rate inside the planet’s Hill sphere.
Methods. We ran global 3D simulations with the gridbased code fargOCA, using a fully radiative equation of state and a dusttogas ratio of 0.01. We built a consistent disk structure starting from vertical thermal equilibrium obtained by including stellar irradiation. We then let a gap open with a sequence of phases, whereby we deepened the potential and increased the resolution in the planet’s neighbourhood. We explored three models. The nominal one features a disk with surface density, ∑, corresponding to the minimum mass solar nebula at the planet’s location (5.2 au), characterised by an α viscosity value of 4 10^{−3} at the planet’s location. The second model has a surface density that is ten times smaller than the nominal one and the same viscosity. In the third model, we also reduced the viscosity value by a factor of 10.
Results. During gap formation, giant planets accrete gas inside the Hill sphere from the local reservoir. Gas is heated by compression and cools according to opacity, density, and temperature values. This process determine the thermal energy budget inside the Hill sphere. In the analysis of our disks, we find that the gas flowing into the Hill sphere is approximately scaled as the product ∑ν, as expected from viscous transport. The accretion rate of the planetary system (envelope plus circumplanetary disk) is instead scaled as √Σv, with its efficiency depending on the thermal energy budget inside the Hill sphere.
Conclusions. Previous studies have shown that pressuresupported or rotationally supported structures are formed around giant planets, depending on the equation of state (EoS) or on the opacity; namely, on the dust content within the Hill sphere. In the case of a fully radiative EoS and a constant dust to gas ratio of 0.01, we find that lowmass and lowviscosity circumstellar disks favour the formation of a rotationally supported circumplanetary disk. Gas accretion leading to the doubling time of the planetary system of > 10^{5} years has only been found in the case of a lowviscosity disk.
Key words: planets and satellites: gaseous planets / protoplanetary disks / planetdisk interactions
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The discovery of PDS70b and PDS70c, the first protoplanets directly imaged in their host disk (Keppler et al. 2018; Haffert et al. 2019), provides a unique opportunity to gain new insights in the processes of planet formation and evolution and how these giant planets acquired their atmospheric composition (Cridland et al. 2023). PDS70b and c are giant planets located in a large dustdepleted cavity, at 21 and 34 au, respectively. Their properties remain unclear but all studies are in agreement with the assumption of high masses (3–12 M_{Jup} Müller et al. 2018; Haffert et al. 2019; Stolker et al. 2020; Wang et al. 2020). Interestingly, both planets are aptly detected in the Hα line and PDS70c is also detected in the UV continuum (e.g. Haffert et al. 2019; Hashimoto et al. 2020; Zhou et al. 2021), indicating ongoing accretion of gas onto the planets and possibly accretion shocks (Aoyama et al. 2018). Spectrophotometric studies show infrared spectral energy distributions, which are notably red and flat, indicating the existence of dust located within the planetary atmosphere or in a surrounding circumplanetary disk (Christiaens et al. 2019; Wang et al. 2021). The presence of small grains around both planets (and in the inner disk within the cavity) is also suggested by hydrodynamical simulations by Bae et al. (2019). Interestingly, evidence for dust grains colocated with PDS70c and b was found with ALMA submillimeter continuum observations (Isella et al. 2019; Benisty et al. 2021). While the emission around PDS70c, closer to the outer protoplanetary disk, is compact and unresolved (at 2 au resolution), the emission around PDS70b, which is located furtherin within the dustdepleted cavity, is faint and diffuse, with a morphology that remains unclear. It is thought that PDS70b is starved of dust grains from the outer disk, as they are filtered out by PDS70c (Pinilla et al. in prep.). While both planets seem to be surrounded by dust, the different observational properties of the emissions around PDS70b and c calls for a better understanding of the structure of the material around the planets and, in general, of the conditions for circumplanetary disks to form.
In this context, numerical simulations can guide our understanding and address major questions such as the formation of a circumplanetary disk (CPD, hereafter), its structure and accretion properties. From the numerical perspective studying the gas dynamics around a giant planet is challenging. As giant planets accrete from the protoplanetary disk and part of the gas enters in the Hill sphere, numerical models require us to resolve both the full disk behaviour and the dynamics in the planet’s Hill sphere. This goal is achieved by using nested meshes on grid based codes either on the full disk or with a local approach, namely, shearing sheet approximation (D’Angelo et al. 2003; Klahr & Kley 2006; Machida et al. 2008; Tanigawa et al. 2012; Szulágyi et al. 2014, 2016), nonregular grid geometry in the planet’s Hill sphere (Lambrechts & Lega 2017; Fung et al. 2019; Lambrechts et al. 2019; Schulik et al. 2019, 2020), or using SPH codes (Ayliffe & Bate 2009).
The first studies have been done in 2D models (D’Angelo et al. 2002) and extended to three dimensions as soon as computing capabilities allowed for it. The reason for including the vertical dimension for giant planets is that their Hill radius is comparable to disk scale height and, in fact, vertical motion has been shown to have an important role for a planet’s accretion (see for example Bate et al. 2003; Klahr & Kley 2006; Tanigawa et al. 2012; Szulágyi et al. 2014; Morbidelli et al. 2014). Comparing the numerical results from the literature is not a simple task, as different codes are used and different physical models are involved. Some of the studies consider the isothermal approximation (Machida et al. 2008, 2010; Tanigawa et al. 2012; Szulágyi et al. 2014; Fung et al. 2019; Béthune & Rafikov 2019), while others adopt a more complex equation of state (EoS hereafter) going from adiabatic (Machida et al. 2010; Fung et al. 2019) to fully radiative EoS on inviscid or viscous disks (D’Angelo et al. 2003; Klahr & Kley 2006; Ayliffe & Bate 2009; Szulágyi et al. 2016; Lambrechts & Lega 2017; Lambrechts et al. 2019; Schulik et al. 2019, 2020; Marleau et al. 2023); in this case, the opacity depends on physical disk properties such as the density and the temperature, or is simply considered as constant (Schulik et al. 2020).
A common result of isothermal models is the formation of a CPD around giant planets. All models agree that the gas that reaches the planetary system (planet plus circumplanetary disk) mainly comes from polar regions and flows down towards the planet with a spiral pattern. A highly debated question among the models is whether gas is also accreted from the equatorial direction. In the case of inviscid disks, some authors have shown that gas reaching the planetary system from disk midplane may lose angular momentum due to the torque exerted by the star through the spiral wave (Szulágyi et al. 2014), whilst others (Klahr & Kley 2006; Machida et al. 2008; Tanigawa et al. 2012) have not found any midplane accretion. However, even in the case of midplane accretion, Szulágyi et al. (2014) found that only 10% of the accreted gas comes from disk midplane, while the remaining accretion was due to polar flow.
In contrast, when considering fully radiative EoS, previous studies are not in agreement: some authors find fully rotationally supported CPD (Ayliffe & Bate 2009; Marleau et al. 2023), whilst others report a slow rotation and a pressure supported planetary envelope (Klahr & Kley 2006; Szulágyi et al. 2016; Lambrechts et al. 2019; Krapp et al. 2024). Using the adiabatic approximation, Fung et al. (2019) showed that a spherical pressure supported adiabatic envelope is found instead of a CPD. Therefore, the isothermal and the adiabatic approximations give the two opposite situations of an extremely efficient and inefficient cooling in the planet’s environment and we can expect to find either a pressuresupported envelope or a CPD when considering the more realistic fully radiative EoS model, depending on the planet’s environment.
Concerning the vertical inflowing gas, it is travelling nearly in free fall at supersonic speed in the cases where a rotationally supported CPD is found (Fung et al. 2019; D’Angelo et al. 2003; Marleau et al. 2023) and subsonic speeds are reported for pressure supported adiabatic envelopes (Fung et al. 2019; Béthune & Rafikov 2019). The inflowing gas hits the envelope’s surface or the CPD and if the shocks occurs at a velocity above a given critical value, the emission of hydrogen H_{α} lines is expected (D’Angelo et al. 2003; Marleau et al. 2023).
The outcome of fully radiative models has been shown to depend on the temperature at the planet’s position (Szulágyi et al. 2016), which can itself be regulated by the opacity (Ayliffe & Bate 2009; Schulik et al. 2019, 2020). More precisely, keeping the gas bound to the planet cold or decreasing the dusttogas ratio to have a low opacity in the local environment enhances the rotational support with respect to the pressure support. The gas accretion rate onto the planet seems independent on the rotational support of CPD: typical accretion rates for a Jupitermass planet take values of the order of 2–6 × 10^{−5} Jupiter masses per year (Bate et al. 2003; Klahr & Kley 2006; Ayliffe & Bate 2009; Machida et al. 2010; Lambrechts et al. 2019; Schulik et al. 2020); with a larger value of 10^{−4} Jupiter masses per year in the case of a low opacity disk found in Schulik et al. (2020). Even the lower value of 2 × 10^{−5} Jupiter masses per year gives a doubling mass time an order of magnitude shorter than typical disk lifetime, enforcing the idea that Jupitermass planets emerged late in disk lifetime. A decline of mass accretion rate is expected for planets embedded in wide gaps (Machida et al. 2010). This could be the case of giant planets embedded in lowviscosity disks.
Recently, Nelson et al. (2023) measured accretion rates on Jupitermass planets by using a simple hydrodynamical model with star accretion on surface layers induced through angular momentum loss. These authors found accretion rates in the interval [10^{−7}, 4 × 10^{−6}] Jupiter masses per year. This result does not require giant planets to form later on. Moreover, the study of lowviscosity disks has garnered interest in recent years both from observational constraints (Pinte et al. 2016; Villenave et al. 2022) and from theories revisiting the generation of turbulence in disks via the magnetorotational instability (Turner et al. 2014). Recent magnetohydrodynamical models (Bai 2017; Béthune et al. 2017) that take into account nonideal effects have shown that in disk with a low bulk viscosity (α ~ 10^{−5}), gas can be transported to the star by thin accretion layers located at the disk surface at rate compatible with the observed staraccretion rates.
The aim of this paper to describe the dynamics and physical structure of the gas in the immediate vicinity of a giant planet and investigate the role of key properties of the protoplanetary disk, such as the mass or the viscosity in determining the outcome of the CPD and the planet’s accretion rate. We follow a selfconsistent approach that goes from the study of the disk structure, whilst also considering the formation of a gap in the protoplanetary disk. This paper is the first of a series of (at least) two papers, with paper II focusing on the chemical delivery of volatile to the planet (Cridland et al. 2024). The paper is organised as follows. In Sec. 2, we present the physics and the numerical models. We also compare global properties of the three disk models with an embedded Jupitermass planet. In Sect. 3, we provide details of the gas dynamics focusing on the Hill sphere. Our conclusions are given in Sec. 4.
2 Models
2.1 Physical model
In this work, the protoplanetary disk (PPD hereafter) is treated as a nonselfgravitating gas whose motion is described by NavierStokes equations. We used the gridbased fargOCA code (Lega et al. 2014) with a specific nonuniform grid setting introduced in Lambrechts & Lega (2017). We used spherical coordinates (r, φ, θ), where r is the radial distance from the star (which is at the origin of the coordinate system), φ is the azimuthal coordinate measured from the xaxis, and θ is the polar angle measured from the zaxis (the colatitude). The midplane of the disk is located at the equator, $\theta =\frac{\pi}{2}$. We worked in a coordinate system that rotates with the angular velocity of a planet of a mass, m_{p}: $${\mathrm{\Omega}}_{p}=\sqrt{\frac{G({M}_{\star}+{m}_{p})}{{r}_{p}{}^{3}}},$$
where M_{*} is the mass of the central star, G is the gravitational constant, and r_{p} is the semimajor axis of a planet assumed to be on a circular orbit. We consider a planet orbiting a Solarmass star and, therefore, we indicate the mass of the star as M_{⊙} in the following.
The dynamics of the gas is provided by the integration of the NavierStokes equations composed by a continuity equation and a set of three equations for the momenta (Lega et al. 2014). To the NavierStokes system of equations, we added an equation for the energy density of the thermal radiation, E_{r}, and that of the internal energy, e = ρc_{υ}T, where ρ and T are the gas volume density and the gas temperature, whilst c_{υ} is the specific heat at constant volume. We followed the evolution of both quantities using the socalled two temperature approach (Commerçon et al. 2011): $$\{\begin{array}{l}\frac{\mathrm{\partial}{E}_{\mathrm{r}}}{\mathrm{\partial}t}\mathrm{\nabla}\cdot \mathit{F}=\rho {\kappa}_{p}({a}_{r}{T}^{4}c{E}_{\mathrm{r}})\\ \frac{\mathrm{\partial}e}{\mathrm{\partial}t}+\mathrm{\nabla}\cdot (e\mathit{v})=P\mathrm{\nabla}\cdot \mathit{v}\rho {\kappa}_{p}({a}_{r}{T}^{4}c{E}_{\mathrm{r}})+{Q}^{+},\end{array}$$(1)
where $\mathit{F}=\frac{c\lambda}{\rho {\kappa}_{r}}\mathrm{\nabla}{E}_{\mathrm{r}}$ is the radiation flux vector calculated in the flux limited diffusion approximation (FLD, Levermore & Pomraning 1981) with λ as the flux limiter (Kley et al. 2009). We indicate with κ_{p} and κ_{r} respectively the Planck and the Rosseland mean opacity, σ is the StefanBoltzmann constant, c is the speed of light. In this paper we consider κ_{p} = κ_{r} (see Bitsch et al. 2013) and use the opacity law of Bell & Lin (1994).
We considered an ideal gas of pressure P with equation of state: P = (γ − 1)e for a gas with adiabatic index of γ = 1.4. The terms P∇ · υ and Q^{+} are the compressional heating and the viscous heating, respectively (Mihalas & Mihalas 1984). We did not include the heating from the central star. This choice allows us to reduce the computational cost with no significant impact on the gas dynamics in the planet vicinity (see Lega et al. 2015). However, the heating from the central star has an impact on the thermal properties of the disk. In Appendix A, we describe the procedure used to find the disk surface temperature that corresponds to the disk aspect ratio in the vicinity if the planet relative to that obtained for a disk heated by the central star.
Disk models parameters.
2.2 Disk models
We modelled a full disk annulus with an azimuthal domain in [−π, π] since we are interested in the analysis of the exchange of gas between the planet’s Hill sphere and the protoplanetary disk. We are also required to trace the gas dynamics on small scales well inside the Hill sphere. For this purpose, we use a nonuniform grid geometry with a grid suitably designed to have small and nearly uniform cells around the planet and large gridcells further out. This procedure was introduced by Fung et al. (2015) and tested for the code fargOCA by (Lambrechts & Lega 2017; Lambrechts et al. 2019). We remark that our small scales inside the Hill sphere correspond to about 10 Jupiter radii (considering Jupiter’s current radius); therefore, we cannot study the phase of growth of Jupiter nor define a boundary between the planet and the surrounding gas. In this paper, we consider three simulations of a disk with constant viscosity and surface density defined as: ∑ = ∑_{0}(r/r_{p})^{−1/2}:
(1) The ‘nominal’ case, hereafter referred to as N, for which we consider $\left({\mathrm{\Sigma}}_{0}/{M}_{\odot}{r}_{p}^{2}\right)=6.67\times {10}^{4}$, corresponding to about 210g/cm^{2} at the planet’s location. For this simulation the viscosity is $\nu /{r}_{p}^{2}{\mathrm{\Omega}}_{p}={10}^{5}$ in code units (the value 10^{−5} in code units corresponds to 10^{11}m^{2}/s or 10^{15}cm^{2}/s and to a value of α = 4 10^{−3} at the planet’s position of 5.2 au). These disk parameters are the same as those reported in Szulágyi et al. (2016); Lambrechts et al. (2019).
(2) The ‘lowmass’ case, referred to as Lm, in which the disk is ten times less massive than the nominal one.
(3) The ‘lowmass, lowviscosity’ case, referred to as LmLv, where we consider a lowmass disk in which we also consider a tentimes smaller viscosity value. In Table 1, we give the simulation names and the main parameters.
2.3 Disks with embedded planets: Numerical procedure
2.3.1 Setup
Before inserting the planet, we ran 2D (r, θ) simulations until the disks reach thermal equilibrium. Thermal equilibrium is obtained with a numerical procedure described in Appendix A. We then expanded the disk azimuthally on φ in the interval [−π, π] and restarted our simulations from initial conditions in thermal equilibrium. We embedded a planet of 20 Earth masses in each disk model and smoothly increased (in a time interval of 20 orbits) the mass until the Jupiter mass was reached. The planet is on the midplane at r_{p} = 5.2au and azimuth φ = 0 (or equivalently at (x_{p}, y_{p}) = (5.2, 0) [au]) on a disk of radial domain [2, 13] [au].
2.3.2 Planetary potential
The gravitational potential of the planets on the disk (Φ_{p}) is modelled as in Kley et al. (2009) using the full gravitational potential for disk elements having distance d from the planet larger than a fraction ϵ of the Hill radius, and a smoothed potential for disk elements with d < r_{sm} ≡ ϵR_{H} according to: $${\mathrm{\Phi}}_{p}=\{\begin{array}{ll}\frac{{m}_{p}G}{d}& d>{r}_{\mathrm{s}\mathrm{m}},\\ \frac{{m}_{p}G}{d}f\left(\frac{d}{{r}_{\mathrm{s}\mathrm{m}}}\right)& d\le {r}_{\mathrm{s}\mathrm{m}},\end{array}$$(2)
with $f\left(\frac{d}{{r}_{\mathrm{s}\mathrm{m}}}\right)=[{\left(\frac{d}{{r}_{\mathrm{s}\mathrm{m}}}\right)}^{4}2{\left(\frac{d}{{r}_{\mathrm{s}\mathrm{m}}}\right)}^{3}+2\frac{d}{{r}_{\mathrm{s}\mathrm{m}}}]$. We recall that the Hill radius is defined as: R_{H} = r_{p}(m_{p}/M_{⊙})^{1/3}. We remark that (unlike the case of 2D simulations) there is no physical need for the use of a smoothed potential in 3D simulations since the full gas column is modelled. However, the planetary potential is singular at the planet’s position and numerically handling the function 1/d when d → 0 requires particular care. Specifically, and as shown in Lambrechts et al. (2019), Fig. A.1, a minimal smoothing length of eight grid cells is required to obtain convergent results with respect to the grid resolution. Here, we adopt the practice suggested by Lambrechts et al. (2019) of doubling the resolution when halving the smoothing length with a multiphase prescription. In the first phase (A hereafter), we use a moderate resolution (N_{r}, N_{φ}, N_{θ}) = (228, 680, 32) on a uniform grid. In this way, we can benefit from the large time step provided by the FARGO algorithm Masset (2000). At this stage, we have about ten gridcells in a Hill diameter and therefore we use a large smoothing length for the gravitational potential: ϵ = 0.8 R_{H}. After the planet’s growth, we continue the simulation for additional 80 orbits to let the gap open. Although the gaps have not converged on this integration time (Schulik et al. 2019), the runaway growth timescale is very short, compared to the viscous timescale required for gaps to reach their steady state. Figure 1 (left column) shows the radial density profiles for the three models at the end of phase A. The density of the Lm and LmLv cases has been multiplied by 10 in order to superpose the background unperturbed profiles and appreciate the differences in the gap’s depth and width. We remark that in the lowmass case, the gap is deeper than in the nominal one because the disk has a smaller aspect ratio (see Fig. A.1 bottom panels). In the low viscosity case, the gap is wider than in the other two models (as expected).
We continue the simulation (phase B) for additional 25 orbits using a regular grid with resolution (N_{r}, N_{φ}, N_{θ}) = (456, 1360, 64) in the three directions. The smoothing length is smoothly reduced with respect to phase A down to 0.4 R_{H}.
2.3.3 Boundary conditions
In the azimuthal direction, we used periodic boundary conditions; whereas in the radial direction, we used evanescent boundaries to minimise the wave reflections (de ValBorro et al. 2006). The reference fields (density, velocity, and gas energy) for evanescent boundaries are those that result from the thermal equilibrium obtained running 2D (r, z) axisymmetric simulations, as described in Appendix A. In the vertical direction at the disk surface, we considered reflecting boundaries for the velocity components. We also copied the density from the neighbours active cells and set the temperature at the values T_{b} (as computed in Appendix A). To limit the computational time, we modelled a single disk hemisphere by using a mirror symmetry at the midplane.
2.3.4 Increase of numerical resolution
In the next phase (C hereafter), we used the nonuniform grid geometry (with a number of gridcells: (N_{r}, N_{φ}, N_{θ}) = (298, 1598, 65)) on a disk annulus of radial domain of r/r_{p} = [0.6, 1.4] (or [3.1, 7.3] [au]) with r_{p} = 5.2 [au]). The grid has finest resolution of $\frac{\mathrm{\Delta}r}{{r}_{p}}=0.0013$ in the three dimensions, close to the planet’s position. The spacing is quasiconstant in the Hill sphere, increases quadratically and has a maximum resolution ratio of 6 in the radial direction (with 10 and 2.5 in the azimuthal and in the vertical directions, respectively, not shown). The number of gridcells in a Hill radius is of 50. We discuss the convergence of results with respect to resolution and smoothing length in Appendix B. The integration annulus fully contains the gap and the radial domain for evanescent boundary conditions, as seen in Fig. 2.
The reference fields for the evanescent boundaries are the gas fields inherited by those computed at the end of phase B. Because of the nonuniform geometry of the grid in azimuth, we cannot use the FARGO algorithm in phase C; thus, the simulations take a longer time and we ran them for an additional ten orbits. The smoothing length in phase C is reduced with respect to phase B down to ϵ = 0.2 R_{H} in 1 orbital period (corresponding to 10 gridcells in the smoothing lenght). All of the figures in the following are snapshots at five orbital periods out of a total of ten orbital periods for this phase. From phase A to phase C, we increased the resolution and reduced the smoothing length: more gas is accreted in the planet vicinity (Fig. 1, top right panels). In Fig. 1 (middleright panel), we observe averaged temperatures similar to the results of phase A on the coarsely grained grid, along with a huge increase (of a factor of 5–20) in the temperature at the planet’s position. This increase in temperature near the planet follows from the higher spacial resolution, along with the deeper potential and the fact that the main source of heating is compression (Lambrechts et al. 2019). The midplane temperature suddenly drops at the boundary of the Hill sphere; however, small temperature spikes are observed near this boundary in the three models due to the intersection of the selected (azimuthal) slice with the spiral wake (as we show in the next sections). We can already see from this study that local properties of temperature and density in the Hill sphere of a giant planet that has opened a gap depend on such disks parameters as the mass and the viscosity.
3 Gas dynamics in the Hill’s sphere
In this section, we provide a detailed description of the dynamics in the planet’s vicinity and quantify the amount of gas that flows through the Hill sphere as well as the accretion rate in the Hill sphere. In a followup paper (Cridland et al. 2024), we will present the chemical evolution of gas along streamlines computed in this work.
To focus on the dynamics in the planet’s vicinity, we worked in the planetocentric reference frame and used either Cartesian coordinates (x, y, z) or spherical coordinates (R, ϕ, ϑ), with R representing the radial distance of the planet, ϕ the azimuthal angle in [0: 2π], where the xaxis, namely, ϕ = 0, marks the Sunplanet direction; ϑ is the colatitude in the interval [0: π/2], with π/2 indicating the planet’s equatorial plane (corresponding to the disk’s midplane) and 0 indicating the pole of the sphere. In Fig. 3, we provide the density field at the disk midplane (left panels) and the azimuthally averaged temperature field (right panels) represented on the plane $(\stackrel{~}{R},z)\equiv (R\mathrm{sin}(\vartheta ),R\mathrm{cos}(\vartheta ))$ for the three models N (top), Lm (middle) and LmLv (bottom). The velocity field (υ_{x}, υ_{y}) is overplotted on the background midplane density with white arrows. On the right panels, the velocity field is azimuthally averaged and density weighted as follows: $$<X{>}_{\rho}=\frac{{\int}_{0}^{2\pi}\rho X\mathrm{d}\varphi}{{\int}_{0}^{2\pi}\rho \mathrm{d}\varphi},$$(3)
where < X > accounts for both components (< υ_{R} >, < υ_{z} >). The colours of the arrows correspond to the vertical component < υ_{z} >_{ρ} in units of the sound speed.
Fig. 1 Radial density profiles (top), azimuthally and vertically averaged temperature (middle), and disk aspect ratio (bottom) profiles for the three models at the end of phase A (left) and C (right). The grey dashed lines indicate the boundary of the Hill sphere and the black dotted line the position of the planet. The density of the Lm and LmLv cases has been multiplied by 10 in order to superpose the density background and appreciate the differences in the gap’s depth and width. We notice that in the LmLv case a large amount of gas is accumulated in a radial domain corresponding to the Hill sphere. This is not gas accumulated around the planet but it is the signature of the libration region around the L4 and L5 Lagrangian points not being depleted yet (a known fact in the case of lowviscosity disks). In the middle panel, the midplane temperature is overplotted with same colours (dashed) on a slice at the planet’s azimuth φ = 0. The higher resolution associated with the deepening of the potential (decreasing smoothing) results in a huge temperature increase up to [1600, 1100, 600] K in the nominal model, the Lm case, and the LmLv model, respectively. In the bottom panel, we remark that it is only in the nominal case that the formation of the gap modifies the aspect ratio considerably. 
3.1 Midplane views
On the plot displaying the density midplane, we clearly see the spiral wake as an overdense region that shrinks and becomes very thin when moving from the nominal to the LmLv case. The central massive region (envelope) extends to about 0.5 R_{H}, gas is subsonic inside the envelope; whereas it approaches the wake at supersonic speed and is then deviated and slowed at subsonic speed after shock with the wake. We remark that the density structure is similar to the one described in Schulik et al. (2020, see their Fig. 3) for similar protoplanetary disk initial conditions (same viscosity and twice the smaller surface density at the planet’s location), but with a small opacity value of 0.01g/cm^{2}. The wake structure seems not to depend on the opacity, however, the region inside 0.4 R_{H} is in our case almost not rotating. This is typical of pressuresupported envelopes and by lowering the opacity, Schulik et al. (2020) found a rotationally supported structure inside 0.4 R_{H} (and strong spiral wakes in the CPD). In the middle and bottom panels (Lm and LmLv cases), in addition to the spiral wake, we observe additional spiral arms: this is a sign that in these cases the gas surrounding the planet forms a rotationally supported structure or CPD with its tidal spiral arms. When the gas approaches the protoplanet, its motion is deviated by the gravity of the planet. From the literature, we have the definition of three regions (Machida et al. 2008, 2010; Tanigawa et al. 2012) corresponding to classical circulation (also called a passby region) and horseshoe motion at distances of >R_{H}, along with a third region for distances of R<R_{H} (called atmospheric region), where the gas is highly perturbed by the planet.
Inside the Hill sphere, we can identify the region where the gas is bound to the planet through energy considerations. Actually, the gas escapes the planet’s potential well if its total energy is larger than the effective potential at the Hill radius: ${\mathrm{\Phi}}_{H}=G{m}_{p}/{R}_{H}3/2{\mathrm{\Omega}}_{p}^{2}{R}_{H}^{2}$; here, the second term is the centrifugal contribution in the local approximation (Hill’s equations). Therefore, we can say that the gas is bounded where E ≤ 0 with the specific energy defined as: $$E=\frac{{v}^{2}}{2}+{c}_{v}T+{\mathrm{\Phi}}_{p}\frac{3}{2}{\mathrm{\Omega}}_{p}^{2}{x}^{2}{\mathrm{\Phi}}_{H}.$$(4)
The zerolevel energy contour (black line in the left panels of Fig. 3) shows that the gas is bounded to the planet up to a distance of about 0.4 R_{H}, with a spherical geometry in the N case and with an elongated or disk like structure on the Lm and LmLv cases.
In the left panels of Fig. 3, we have also plotted forward (full lines) and backward (dashed lines) integration for the set of initial conditions of Table 2 projected on the midplane (z = 0), using the 2D velocity fields components (υ_{R}, υ_{&phis}). Ignoring the vertical motion, the flow pattern emphasises the influence of the wake on the gas motion. In the top panel (nominal case), after being deviated by the shock with the spiral wake, the brown and pink lines according to the velocity field are further deflected near the Lagrangian point, L_{1}, at approximately x = −0.6, y = 0.1 R_{H}, and reach the boundary of the hot envelope at about 0.4 R_{H}. On average the radial motion of the midplane is directed inward, gas moves towards the hot central envelope and, in the full 3D case (Fig. 3, topright) is then lifted up. A consequent depression in the midplane will keep the gas moving toward the hot envelope (“green” and “magenta” streamlines), while the others (‘cyan’, ‘black’, and ‘orange’) pass further away from the planet and have a typical horseshoe motion.
In the Lm case, the brown and pink lines are deviated by the wake, have about half a rotation around the planet and are bent again by respectively the CPD spiral arm and the wake and finally leave the Hill region, while the green and magenta lines spiral out and are kinked at each rotation by both the wake and the CPD spirals. A similar behaviour is observed in the LmLv case, with the majority of streamlines entering the Hill sphere participating to horseshoe motion.
Fig. 2 Vertically integrated density of the simulated nominal disk (simulation N in Table 1) in the highresolution phase (phase C, see text) with nonuniform grid geometry. We clearly see the spiral arms originating at the planet position (x_{p}, y_{p}) = (5.2, 0) [au], revealing the interaction between the planet and the disk. The disk annulus fully contains the gap and slightly extends over the edges. 
3.2 Vertical views
The contour levels plotted in Fig. 3, right panels correspond to the vertical component of the specific angular momentum, j, in a non rotating frame centered on the planet, normalised over the Keplerian angular momentum relative to the planet. More precisely, we have: $$j=\frac{<{v}_{\varphi}\stackrel{~}{R}{>}_{\rho}+{\mathrm{\Omega}}_{p}{\stackrel{~}{R}}^{2}}{\sqrt{G{m}_{p}\stackrel{~}{R}}},$$(5)
with $\stackrel{~}{R}$ being the cylindrical radius. In the nominal model, along the polar directions, the gas has less than 30% of the Keplerian angular momentum and flows quasivertically towards the planet with a subsonic sound speed – as previously reported in the adiabatic approximation (Fung et al. 2019; Béthune & Rafikov 2019). The infall speed further decreases when the gas approaches the hot envelope at about z ~ 0.5[R_{H}]. The gas is bound to the planet inside a spherical envelope of radius 0.4 Hill radii (black line showing the zero energy contour, according to Eq. (4) in Fig. 3, topright panel). This is typical of a pressuresupported envelope in agreement with the finding of Szulágyi et al. (2016); Lambrechts et al. (2019); Fung et al. (2019). Within the envelope, the gas speed is negligible with respect to the values observed at larger distances.
When looking at the Lm, LmLv cases, we see that gas flows quasivertically at supersonic speed in the region bound by the line with 0.3 Keplerian angular momentum in the Lm case and on a larger fraction of the Hill sphere in the LmLv case. In the Lm case, we observe a dense quasispherical envelope of a radius of ~0.2[R_{H}] and a thick disk that is partially supported by pressure. The separation between the envelope and the surrounding disk can approximately be set at ~0.2[R_{H}]: for larger distances, the angular momentum gets larger than 60% of the Keplerian value. In the LmLv case, the rotational support increases to 80% of the Keplerian one at a distance of $\stackrel{~}{R}$ > 0.2[R_{H}]. In the region bounded to the planet and delimited by the zero energy contour (black line in the right panels of Fig. 3), the averaged motion has negligible velocity with respect to the rest of the Hill sphere. The disk structure appears to shield the gas entering in the Hill sphere at low altitude by deviating it, lifting it up (middle right panel of Fig. 3) above the disk structure, and then moving it radially towards the planet, where it falls almost vertically when the angular momentum is low. We describe the dynamics in the planetary system in Sec. 3.4.
Fig. 3 Parameter details for all three models. Left: midplane density field interpolated in cartesian planetocentric coordinates(x, y) = (R cos(ϕ), R sin(ϕ)) in units of the Hill radius). The cartesian components (υ_{x}, υ_{y}) of the midplane velocity field are plotted with arrows. The coloured lines are straightforward (full) and backward (dashed) integration of streamlines obtained with the 2D midplane velocity field. Right: azimuthally averaged density field interpolated in planetocentric coordinates and represented in $(\stackrel{~}{R},z)=(R\mathrm{sin}(\vartheta ),R\mathrm{cos}(\vartheta ))$ plane. The azimuthally averaged and density weighted velocity is represented with arrows. The colours of the arrows correspond to values of the of the gas speed normalised over the sound speed(respectively the midplane speed (left column) and the absolute value of υ_{z} (right column)). The yellow contour levels correspond to the vertical component of the angular momentum in units of the Keplerian momentum with respect to the planet. The coloured markers give the initial conditions for the streamlines of Table 2, projected on the (x, y) and ($\stackrel{~}{R}$, z) planes, respectively, and the black solid lines identify the 0 energy contour (see Eq. (4)). 
3.3 Rotational support
The fraction of pressure support with respect to rotational support of a planetary system depends on the cooling efficiency in the Hill sphere (Fung et al. 2019). Ayliffe & Bate (2009) had previously shown that the rotational support increases when decreasing the opacity and this result was recently confirmed in a detailed study by Schulik et al. (2020). Here, we use the Bell & Lin opacity k(ρ, T) for our three models with a dusttogas ratio of 0.01. Although the inner and denser part of the planetary envelope may be convective (Zhu et al. 2021), in the outer part, or in the CPD forming region, cooling is mainly radiative (Fung et al. 2019; Krapp et al. 2024). Radiative cooling can be quantified by computing the timescale over which radiation is diffused over the Hill sphere. Following the definition of Bitsch & Kley (2010) we compute: $$\tau =\frac{{R}_{H}^{2}}{{D}_{r}},$$(6)
where D_{r} is a radiative diffusion coefficient in the optically thick limit: $${D}_{r}=\frac{4ca{T}^{3}}{3{c}_{v}{\rho}^{2}k(\rho ,T)}.$$
In Fig. 4 we show, as a function of the radial distance from the planet, quantities averaged on spherical planetocentric shells, according to: $$F(R)=\frac{1}{2\pi {R}^{2}}{\int}_{0}^{2\pi}{\int}_{0}^{\pi /2}f(R,\theta ,\varphi /2){R}^{2}\mathrm{sin}(\vartheta )\mathrm{d}\vartheta \mathrm{d}\varphi ,$$(7)
where the function, f, is the cooling timescale in units of orbital periods τ/T (topleft), the opacity (top right), and the volume density and temperature (bottom panels).
The inner region up to 0.1 R_{H} has a high density and opacity values larger than 1 cm^{2}/g for the three models (see Fig. 4) with a moderate decrease of the opacity inside 0.1 R_{H}, where the temperature is higher than 1000° and the opacity decreases because of gas ionisation processes. The radiative cooling timescale is of the order of 10^{3} orbital periods for the three models. At distances greater than 0.3 R_{H} from the planet, τ decreases to values smaller that the dynamical timescale for both the Lm and LmLv cases. In a recent paper, Krapp et al. (2024) proposed a necessary condition for the formation of a CPD depending on the value of the radiative cooling timescale. According to this condition, a CPD forms when the cooling time is of the same order or smaller than the orbital period at the tidal truncation radius: ≃0.3 R_{H}. Although we do not have the same definition of τ as in Krapp et al. (2024), this condition is fullfilled by the LmLv case and marginally fulfilled in the Lm case where the disk structure is thick and partially supported by pressure.
3.4 Streamlines
A detailed description of gas dynamics in the planet’s vicinity is given via the description of streamlines in Machida et al. (2008); Bitsch & Kley (2010); Tanigawa et al. (2012) in the case of isothermal and adiabatic models. With the aim of comparing the gas dynamics in our three models. we consider in this section a set of initial conditions on the Hill sphere (see Table 2 and Fig. 3 for projections on the disk midplane and on the vertical slice).
Streamlines give the motion of gas parcels for quasistationary velocity fields. We checked this requirement by comparing the velocity fields at different snapshots on our hydrodynamical integration, which cover nine orbits of Jupiter around the star. For the streamlines, we integrated the velocity field up to ten Jupiter orbits forwards and to two Jupiter orbits backwards. Figure 5 reports the streamlines on a portion of the space slightly larger than the Hill sphere for an integration time of about four Jupiter orbits. On the Hill sphere, we have plotted the radial gas flow rate in planetocentric coordinates: negative values correspond to gas inflow and positive values correspond to outflow (red and blue points, respectively, in Fig. 5). From a simple consideration of the flow sign, we see that gas inflow dominates quasi totally in the Lm and LmLv cases (red regions), whereas a large fraction of the sphere is dominated by outflow (blue regions) in the nominal case. We provide quantitative values of inflow and outflow rate on the Hill sphere as a function of z/R_{H} in Sec. 3.5.
As a first global comment, we can observe that in the N case gas entering at high latitude tends to remain trapped and rotate around the planet, either exiting at high altitude or slowly approaching the planet’s envelope. In the other two cases, a larger number of streamlines are trapped in the planetary system and a more flattened disk structure is observed for the LmLv case, with respect to Lm one. We remark that a larger internal energy budget inside the Hill sphere as for the N case (see Fig. 4 bottom panels) seems to have a key role in determining the fraction of gas accreted in the planetary region with respect to the total gas entering in the Hill sphere. We provide a quantitative study in Sec. 3.5.
In the following, we describe the behaviour of single streamlines by also looking at quantities reported on Fig. 6, where we can appreciate the projected path of streamlines on the x − y global disk as well as the evolution with time of their z and radial components relative to the planet. The dashed lines on the left panels of Fig. 6 correspond to backward integration: we see that in all cases, the gas approaching the Hill sphere comes from Uturn motion passing close to the planet; namely, close to Lagrangian points L1 and L2, from the inner or outer branch of the horseshoe region (in agreement with the findings of Machida et al. 2008, 2010).
Fig. 4 Parameter details for all three models. Topleft: radiative diffusion timescale (from Eq. (6)). Topright: opacity. Bottomleft: volume density. Bottom right: temperature. All the quantities are obtained through averages over spherical shells centred on the planet. 
3.4.1 Nominal case
In the nominal case, the initial conditions that are very close to the planet’s pole (‘dark blue’ and ‘green’ streamlines) gas enter deep into the envelope and is accreted (see also Fig. 6, topmiddle panel), although very low density is carried by gas at high z values. We consider that the gas carried by a streamline is accreted when it enters the region gravitationally bound to the planet (shown by the black lines in Fig. 3). Not all of the gas entering the Hill sphere close to the planet’s pole is accreted, in fact we can see that the ‘magenta’ line, which has an initial value of z/r_{H} very close to the ‘green’ one, but different x and y positions (see Fig. 3, left panel) escapes from the Hill sphere. Precisely, in this case, the gas has a short residence in the Hill sphere moving at lower z values and then is lifted higher up and leaves the Hill sphere reaching the outer disk region (Fig. 6, top row). The computation stops when a streamline reaches the computational domain’s boundaries. Starting further down at z/r_{H} = 0.6 and z/r_{H} = 0.35 (respectively, ‘cyan’ and ‘black’ streamlines), the gas is lifted up in z, according to the direction of the vectors in Fig. 3 (topright panel). In the longer integration presented in Fig. 6, top panels, the gas transits to the outer disk (‘cyan’ streamline), whereas it reenters the Hill sphere (‘black’ streamline) after a horseshoe orbit and, finally, transits into the inner disk (Fig. 6, top left panel). A similar behaviour occurs with respect to the ‘orange’ streamline starting at z/r_{H} = 0.1: the gas makes a Uturn in the Hill region (similar to the midplane streamline of Fig. 3, topleft) followed by an horseshoe orbit, then again approaches the Hill sphere, entering at higher z values at about t = 8 orbital periods; finally, it leaves the horseshoe region and reaches the disk’s upper boundary (see Fig. 6, top panels).
For streamlines with initial conditions that are very close to the midplane, we observe: a close approach to the bounded region at d ≃ 0.4R_{H} and a few rotations around the planet that are similar to the corresponding midplane streamline (‘pink’ line). Then the gas is lifted up and transits to the outer disk remaining at large z values (Fig. 6, top panels). Gas starting on the midplane, but at a x y position close to the L1 Lagrangian point (‘brown’ streamline), is deflected and remains bounded at a distance of about 0.4 Hill radii; namely, at the planet’s envelope boundary (Fig. 6, topright panel). We notice that the dynamical behaviour is similar to the one observed for the corresponding midplane streamline of Fig. 3 (top left); however, in the full 3D velocity field, the gas is lifted up to about 0.4 Hill radii (Fig. 6, top middle panel) and then oscillates vertically whilst rotating around the planet (as indicated by the small period oscillations in Fig. 6, top (middle and right panels). If the gas entering in this case from the midplane carries small dust particles, it will then pollute the planet’s envelope.
Fig. 5 3D representation of streamlines with initial conditions on the Hill sphere (see Table 2). The colour bar corresponds to planetocentric mass flow rate on the Hill surface. The sign is negative for mass inflow and positive for mass outflow. The total integration time is short (about four Jupiter orbits) in order to make this plot clear. The evolution of the same streamlines on a longer integration allows us to follow the dynamics on a projected plane (x/R_{H}, y/R_{H}) and to appreciate the evolution of the distance with respect to the planet in Fig. 6. 
3.4.2 The Lm and LmLv cases
For the other two cases, the dynamics appear very different: the majority of streamlines are trapped in a toroidal disk structure with an outer bound at about 0.4 R_{H} (slightly larger than the tidal truncation radius) and an inner bound at about 0.25 R_{H}. To explain the path of such streamlines, we show the velocity field for distances <0.5R_{H} in the ($\stackrel{~}{R}$, z) plane with normalised vectors (Fig. 7). The vector colours correspond to the norm of $({v}_{\stackrel{~}{R}},{v}_{z})$. The norm of the vectors cover a range of about five orders of magnitude, with smaller values inside the planet’s envelope where the gas is bounded. Gas that shocks the envelope (‘dark blue’, ‘green’, and ‘magenta’ streamlines) at about z ~ 0.2 and z ~ 0.12 (for the (Lm and LmLv cases, respectively) is deviated and trapped in circulating motion in the ($\stackrel{~}{R}$, z) plane (toroidal motion in the full 3D field). Initial conditions making a Uturn far from the planet as the ‘yellow’ and ‘black’ streamlines undergo rapid transit, as observed for the corresponding midplane streamlines of Fig. 3 (middle and bottomleft panels). The ‘brown’ streamline is deflected by the wake and enters in the inner protoplanetary disk in the LmLv case, whereas it is captured in the disk structure in Lm case (Figs. 5 and 6, middle and bottom panels). When we observe the ‘brown’ streamline, as computed with the midplane velocity field in the Lm case, we see two important deviations when the gas encounters the dense spiral wake, which eventually makes it leave the Hill sphere. Using the full 3D velocity field, the encounter with the wake also modifies the motion in the vertical direction and lifts up the gas (as described in Schulik et al. 2019, Fig. 14, and explained as baroclinic kicks). At the encounter with the outer wake, the streamlines are deviated again and appear to be finally trapped in the disk region. A similar behaviour is observed for the ‘pink’ streamline in the Lm case. Despite the complexity of 3D motion in the Hill sphere, we can conclude by saying that the accretion of gas entering in the Hill sphere is favoured in the Lm and LmLv cases. Here, the gas is colder outside the planet’s envelope with respect to the nominal case, which has a larger internal energy budget inside the Hill sphere (see also Fig. 4, bottom panels).
Fig. 6 Projection of the streamlines of Fig. 5 on the planetocentric x − y plane (left). The blue contours give the boundaries of the disk. Gas approaching the Hill sphere can have a rapid transit phase in the horseshoe region and be released in the inner or the outer disk trapped in the planetary system. Vertical evolution of streamline with respect to time (middle). Evolution of the distance with respect to the planet as a function of time (right). The dashed lines correspond to backward integration: in the left panel we see that gas approaching the Hill sphere comes from Uturn motion passing close to the planet from the inner or outer branch of the horseshoe region. 
Fig. 7 Same as Fig. 3, right panels showing a zoom close to the planetary system. Vectors of components $({v}_{\stackrel{~}{R}},{v}_{z})$ are normalised and the colours correspond to the norm values. The norms cover a range of about five orders of magnitude, with smaller values deep inside the planet’s envelope. 
Fig. 8 Absolute value of the inflow and outflow mass rate as a function of z on the Hill sphere. In the Lm and LmLv cases the mass inflow is larger than the mass outflow except very close to the midplane which appears decreating. Instead, in the N case, the Hill sphere is accreting up to z = 0.2[R_{H}] and for z > 0.7[R_{H}]. In all the models, the contribution to the total mass flow in the Hill sphere becomes negligible at z > 0.8[R_{H}]. 
Total mass inflow, $\left{\dot{M}}_{i}\right$ and outflow, ${\dot{M}}_{o}$ rates.
3.5 Accretion
From the values of mass flow rate (plotted on the Hill sphere in Fig. 5), we have a qualitative view of regions of gas inflow and outflow. We sum these values on the azimuthal direction to provide the total inflow ${\dot{M}}_{i}(z)$ and outflow mass flow rate ${\dot{M}}_{o}(z)$ as a function of z[R_{H}].
The resulting values are plotted in Fig. 8, taking the absolute value for the inflow rate. The difference between the inflow and the outflow mass rate will give the net accretion rate, simply called the ‘accretion rate’ in the following. In the Lm and LmLv cases (Fig. 8), the mass inflow is larger than the mass outflow except very close to the midplane. In other words, very close to the midplane the flow is decreating and the rest of the Hill sphere is accreting. Instead, in the N case, the Hill sphere is accreting up to z = 0.2[R_{H}] and for z > 0.7[R_{H}]. In all the models, the contribution to the total mass flow in the Hill sphere becomes negligible for z > 0.8[R_{H}]. The total mass flow in and out of the Hill sphere and the resulting accretion rate are reported in [M_{⊕}/yr] in Table 3.
We can now quantify the accretion rate inside the Hill sphere using the whole time series of data at our disposal. More precisely, we first compute the total mass m_{R}(t) contained on semispheres of radius, R, centred on the planet with R going from the 0.01 to 1 Hill radius. We have a total integration time of ten orbits (phase C; see Sec. 2) and the first orbit is used to smoothly reduce the smoothing length. Therefore, we obtained the mass accretion rate from a linear interpolation of the function 2m_{R}(t) (where the factor of 2 gives the whole mass on the sphere) over the remaining nine orbits. In Fig. 9, we plot the accretion rate in Earth masses per year. At the Hill radius, we also plot the total rate of inflow and the accretion rate from Table 3. We observe the good agreement among the accretion rates computed through radial flow on a single time snapshot and the one obtained through interpolation over nine orbits.
In Fig. 9, we notice that below 0.3 Hill radii for the nominal case (0.2 and 0.1 Hill radii for the Lm and LmLv cases, respectively), the accretion rate suddenly decreases and becomes negligible when approaching the planet. We can consider the interior envelope to be slowly contracting and accreting on larger timescales, with respect to our short integration time.
The accretion rate reaches a constant maximal value close to the boundary of the zero energy region (see Fig. 3 and Eq. (4)). This means that gas accreted by the planet piles up in the outer envelope in the N case, and in the CPD in the other two cases; in agreement with our dynamical study.
We observe that in the nominal case, the rate of inflow is of ~0.1M_{⊕}/yr (Table 3, in agreement with Lambrechts et al. (2019, Fig. 2) whilst the accretion rate is of ~0.029M_{⊕}/yr. This is due to the large exchange of gas in the Hill sphere and corresponds to the qualitative description we described in the previous section, via the streamlines.
Instead, in the Lm case, the rate of inflow (~0.011M_{⊕}/yr, Table 3) is marginally higher than the accretion rate (0.009M_{⊕}/yr, Table 3). Therefore, although there is a difference of about a factor of 10 in the total inflow rate between the nominal and the Lm case, as we expect for a viscous transport proportional to ∑_{ν}, there is only a factor of 3 in difference for the accretion rate between the two cases. This result was not expected, since the accretion is considered to behave proportionally to the mass inflow rate; however, we have seen that the accretion depends on the dynamics in the Hill sphere and this changes drastically from the nominal case in comparison to the other two cases, which appear to accrete the in flowing gas more efficiently. Precisely, the accretion appears to be efficient also in the LmLv case, with an accretion rate of 0.0015M_{⊕}/yr for a total inflow rate of 0.0023M_{⊕}/yr. From the viscous disk evolution, the quantity of gas that may be supplied by the disk and can flow onto the planet is proportional to ∑_{ν}. The mass inflow on the Hill sphere provided in Table 3 actually scales as expected based on viscous evolution. Instead, for the net accretion rate in the Hill sphere, we have obtained a scaling $\propto \sqrt{\mathrm{\Sigma}}\nu $.
The accretion rate inside the Hill sphere is often used to infer the mass doubling time of the planet (or of the planetary system when considering the envelope and the circumplanetary disk). In our case, the Jupitermassdoubling time is of the order of a few 10^{4} of years, as previously reported in the literature, for both the N case and the lowmass disk case. Therefore, the typical scenario of a Jupitermass planet entering runaway accretion late in the lifetime of a protoplanetary disk remains valid when considering disks less massive than MMSN at the planet’s location.
In the case of a lowviscosity disk, the mass doubling time increases by an order of magnitude. Thus, it is only in the lowviscosity case that a giant planet may enter in runaway accretion earlier in the disk’s lifetime to achieve the same final mass.
Fig. 9 Accretion rate as a function of the distance to the planet in the Hill sphere computed interpolating the time series of the mass contained in a sphere of radius R over 9 planet’s orbital periods (see text). For each case, we report the rate of inflow (circles) and the accretion rate on the Hill sphere (triangle), as computed on a snapshot at five orbital periods (see Table 3). 
4 Conclusion
Using 3D hydrodynamical simulations with a fully radiative equation of state, we have studied the dynamical properties of gas in the Hill sphere of a Jupitermass planet in three different protoplanetary disks configurations. The first features a disk with a surface density corresponding to the minimum mass solar nebula (at the planet’s location), characterised by constant viscosity, equivalent to having an alpha parameter of 4 × 10^{−3} at the location of the planet. The second disk model differs from the nominal one in mass, with a surface density that is ten times smaller than the nominal one and same viscosity; whilst in the third model, we also reduce the viscosity value by a factor of 10.
During gap formation, giant planets accrete gas inside the Hill sphere from the local reservoir. Gas is heated by compression and cools according to opacity, density, and temperature values. This process determines the thermal energy budget inside the Hill sphere. We can then describe the characteristics of the planetary system and the dynamics in the Hill sphere when the gas has reached a quasistationary state.
In the nominal case, a pressuresupported envelope forms, as already reported in the literature (Ayliffe & Bate 2009; Szulágyi et al. 2016; Fung et al. 2019) and gas cools on timescales that are larger than the orbital period on the whole Hill sphere. In the other cases, we observe an envelope and a CPDlike structure supported partially by pressure and partially by rotation. A necessary condition for the formation of a rotational supported structure was recently given in Krapp et al. (2024) based upon the requirement that the cooling time scale at the tidal truncation radius must be comparable to (or shorter than) one orbital period at the planet’s location. We show that our disk model with low surface density and and low viscosity satisfies this criterion, whilst the low surface density case marginally satisfies it; here, the disk structure in the region 0.2–0.4 R_{H} is indeed quite heavily supported by pressure. The dynamics inside the Hill sphere is also shown to depend on the total energy budget, with the gas being accreted preferentially from midplane in the nominal case and lifted up to colder layers, leaving the Hill sphere preferentially at an altitude of 0.3 < z < 0.7 R_{H}. In the other cases, the gas enters in the Hill sphere and accretes on the planetary envelope except for z values very close to the midplane where the radial motion of the gas (relative to the planet) is positive. The polar accretion channel exists in all the cases although a very small fraction (5% in the N case and 1% in the Lm and LmLv cases) of the accreted gas comes from polar directions with z > 0.8 R_{H}.
In all the cases, a hot envelope is formed in the planet’s vicinity extending to ~0.4 R_{H}, ~0.2 R_{H} and ~0.1 R_{H}, respectively, in the (1) nominal, (2) lowmass, and (3) lowmass+lowviscosity cases. A similar result of giant planets forming an envelope (and eventually a disk structure) that is limited to a fraction of the Hill sphere was also found in Schulik et al. (2020); Krapp et al. (2024).
To form a rationally supported structure or CPD for our nominal disk parameters it is necessary to reduce the opacity (Schulik et al. 2020). In a reduced opacity environment the building blocks for satellite formation may be planetesimals captured in the CPD (Ronnet & Johansen 2020). In this paper, we kept the dust to gas ratio to usual circumstellar disk values and considering that, although dust is filtered at the planet’s outer gap edge, particles coupled to the gas are free to cross the gap; thus, this may allow the opacity to remain at disk values (Morbidelli et al. 2023).
Although we are limited in resolution and a finite smoothing length, we observe that by increasing the resolution, the inner envelope becomes hotter and denser; at the same time, the properties in the outer part of the Hill sphere remain basically unchanged, so that our results seem to be robust with respect to resolution and smoothing length (Appendix B and Szulágyi et al. 2016). Concerning the accretion rate and the planet’s mass doubling time, our results are consistent with the literature in the nominal and low mass cases.
In particular, we have shown that the scaling of the net accretion rate among the three models behaves as $\sqrt{\mathrm{\Sigma}}\nu $. Therefore, only reducing the mass does not significantly increase the massdoubling time. Lowviscosity disks instead favour a significant increase in the massdoubling time.
When considering the lowviscosity disk case, we must remark that viscous transport gives a star accretion rate of ~10^{−9} M_{⊙} per year, which is low with respect to observed values of 10^{−7}, 10^{−8} M_{⊙} per year (Hartmann et al. 1998; Manara et al. 2016). Star accretion on surface layers, as provided by MHD models, offers a mechanism for disks with low bulk viscosity to have star accretion rates consistent with the observed values. Recently, Nelson et al. (2023) measured accretion rates on Jupitermass planets by using a simple hydrodynamical model with star accretion on surface layers induced through angular momentum loss (with a rate of ~10^{−8} M_{⊙} per year in magnetised winds). These authors monitored the gas entering in the Hill sphere of a Jupitermass planet by varying the thickness of the accretion layer and found a massdoubling timescale going from 0.2 to 10 Myrs with larger values for very thin accretion layers. The value obtained in this work for the Jupiter doubling time of the planetary system in the LmLv case is of about 0.2 Myr, consistent with the case of gas radially transported across the whole column density in the (Nelson et al. 2023). Applying the procedure presented here on models with thin accretion layers is beyond the scope of this paper and will be undertaken in a future study.
In a followup paper (Cridland et al. 2024) we explore the chemical evolution of the gas along streamlines and the resulting chemical structure of the gas in the vicinity of the planet. The aforementioned dependence of the gas temperature on the spatial resolution could change the chemical structure of the gas in the planet’s neighbour and modify the energy content of any molecular line emission that would emanate from near the planet. However, when reaching temperatures above 2000 K by increasing the resolution and further deepening the potential, the model should take into account hydrogen dissociation, which is not included in this work. One example of an equationofstate adapted to take into account the endothermic reaction of hydrogen dissociation was done through the reduction of the adiabatic index by Fujii et al. (2019). The authors observed temperatures exceeding 10^{4} K in the few grid cells with very deep gravitational potential in cases where a CPD is formed. This is an indication that performing higher resolution studies with a deeper potential would mainly affect the innermost regions close to the planet. However, we consider that further numerical experiments with improved models that take into account hydrogen dissociation, as well as observations, are needed. This will help to better constrain the properties of the gas surrounding embedded planets.
Acknowledgements
We thank an anonymous referee for constructive suggestions that helped improving the manuscript. M.B. and A.C. received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement No. 101002188). A.M. is grateful for support from the ERC advanced grant HolyEarth N. 101019380. This work was granted access to the HPC resources of IDRIS under the allocation 2023 A0140407233 and 2024 A0160407233 made by GENCI. E.L. wish to thank Alain Miniussi for maintenance and refactorisation of the code FARGOCA.
Appendix A Thermal structure
Disks that differ in mass and viscosity have also a different vertical thermal structure. To obtain a coherent description of the thermal structure, we included stellar irradiation in our models using the numerical prescription provided in Bitsch et al. (2014) and performed 2D (r, z) axisymmetric simulations. We considered a solarmass star with a surface temperature of 4370°K, a radius of 1.5 solar radii, and a disk annulus that extends in the interval r/r_{p} = [0.4, 9] (or r = [2.08, 45.9] au). In the vertical direction the disk extends by 20° over the midplane. We considered a resolution of (N_{r}, N_{θ}) = (528, 82).
In our 3D simulations we did not consider stellar irradiation (Eq. 1), in order to have reasonable computational times and to limit the amount of CPUs. However, we can mimic the thermal equilibrium structure obtained in disks heated by the star, by suitably fixing the disk surface temperature.
To this purpose, we ran test cases in which we imposed different values of the boundary temperature, T_{b}, at the disk surface. The simulations have a resolution of (N_{r}, N_{θ}) = (228, 32) on a disk annulus extending radially in the interval r/r_{p} = [0.4, 2.5] and vertically by about 7° over the midplane^{1}.
Figure A.1, left shows the equilibrium midplane temperature (top) and the disk aspect ratios (bottom) for the nominal N model including the stellar heating and for the four test cases without stellar heating. In the inner part of the disk, the values of the aspect ratio are very similar for disks with and without stellar heating (see also Bitsch et al. 2013). This means that viscous heating dominates over the stellar irradiation up to 7–8 au. Beyond 10 au, the irradiated disk is flared, while the disk with T_{b} = 10 °K collapses to small H/r at large radii, as is typical of disks without stellar heating. We mimicked the flaring by imposing a temperature at the disk’s surface higher than T_{b} = 10 °K. For the nominal model, a value of T_{b} = 40 °K nicely fits the temperature structure of the stellar irradiated disk. We remark that the planet position in the paper (vertical black dotted line in Fig. A.1) is at the boundary of the region dominated by viscous heating.
Fig. A.1 Midplane temperature (top) and disks’ aspect ratio (bottom) as a function of radial distance from the host star for 2D simulations for the nominal and Lm cases (left and right, respectively). The red lines include the effect of stellar irradiation. The other simulations have no stellar heating and differ in the boundary temperature, T_{b}, at disk’s surface. The value of T_{b} that allows us to recover the profile of the stellar irradiated disk is used for 3D simulations of in the rest of the paper. Horizontal dotted lines in the top plots correspond to the values of T_{b} in the simulations. The vertical black dotted line marks the position of a Jupitermass planet in the paper. 
When considering the Lm disk (Fig. A.1, right), the viscous heating, proportional to ∑, dominates only in a very narrow region at the inner disk boundary (Fig. A.1, bottom right), and the planet’s location is in the part of the disk where stellar heating dominates. The boundary temperature required to mimic stellar irradiation is T_{b} = 30 °K. In the model with low viscosity (LmLv), the stellar heating dominates in the full disk range that we are considering. We remark that the amount of energy absorbed in the disk from the star’s photons depends on the density and the opacity and therefore the Lm and LmLv disks have the same thermal structure^{2}. Therefore, we also considered a value of T_{b} of 30 °K in the lowviscosity case.
Appendix B Dependence of results on resolution and smoothing length
This appendix presents a study of the dynamic of the gas in the planet’s vicinity for the low mass case by increasing the resolution and by decreasing the smoothing length in order to check the convergence of results with respect to both quantities: the resolution and the smoothing length. For this, we further refined our grid on the same domain previously studied (phase C) and restarted the code from the output at five orbital periods of phase C for two additional orbits on a grid of (N_{r}, N_{ϕ} N_{θ}) = (498, 1998, 60) gridcells for two cases, shown in Fig. B.1:
We kept a smoothing length of ϵ = 0.2. Results (red curves in Fig. B.1) are shown overplotted with respect to those of the Lm case (presented in this paper) at the resolution of phase C (green curves) and showing convergence with respect to resolution.
The new grid has ten grid cells inside 0.1 Hill radii. Therefore, according to the criterion from Lambrechts et al. (2019), by requiring eight grid cells in the smoothing length we can set ϵ to 0.1 R_{H}. On the first restarted orbit we smoothly reduced the smoothing length and we show results (blue curves in Fig. B.1) at the end of the second orbit when quasistationary fields are obtained.
The nominal case was analysed in Szulágyi et al. (2016) using nested meshed with the Jupiter’s code and by simulating disks sectors with the code fargOCA. The resulting gas structure around the planet showed a spherical envelope, like the one that we obtain here at lower resolution, with temperatures increasing in the innermost gridcells when increasing the resolution and decreasing the smoothing length, leaving almost unchanged the density and temperature in the outer part of the Hill sphere. The same conclusion can be drawn here: we notice that the temperature increase inside 0.1 R_{H} with a corresponding decrease in the opacity, which results in a decrease in the radiative timescale. The planetary envelope results to be denser but the profiles for distances larger than 0.1 R_{H} remain unchanged. Therefore, we conclude that our results are robust with respect to both the resolution and the smoothing length.
Fig. B.1 Same as Fig. 4 but for the Lm case (green curves) and for higher resolution simulations with: a smoothing length ϵ = 0.2 (red curves) and ϵ = 0.1 (blue curves). We remark that decreasing the smoothing length impacts the innermost regions (x < 0.1R_{H}) leaving unchanged the outer part of the Hill sphere. 
References
 Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657 [Google Scholar]
 Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJ, 884, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N. 2017, ApJ, 845, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
 Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJ, 916, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Béthune, W., & Rafikov, R. R. 2019, MNRAS, 488, 2365 [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., & Kley, W. 2010, A&A, 523, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & DobbsDixon, I. 2013, A&A, 549, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Morbidelli, E., L., & Crida, A. 2014, A&A, 564, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJ, 877, L33 [Google Scholar]
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cridland, A. J., Facchini, S., van Dishoeck, E. F., & Benisty, M. 2023, A&A, 674, A211 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cridland, A. J., Benisty, M., & Lega, E. 2024, A&A, submitted [Google Scholar]
 D’Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647 [CrossRef] [EDP Sciences] [Google Scholar]
 D’Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [Google Scholar]
 Fujii, Y., Gressel, O., Tomida, K., & Ziegler, U. 2019, EPSCDPS Joint Meeting 2019, 2019, 1408 [Google Scholar]
 Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101 [Google Scholar]
 Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
 Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [Google Scholar]
 Hashimoto, J., Aoyama, Y., Konishi, M., et al. 2020, AJ, 159, 222 [Google Scholar]
 Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25 [Google Scholar]
 Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klahr, H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krapp, L., Kratter, K. M., Youdin, A. N., et al. 2024, arXiv eprints [arXiv:2402.14638] [Google Scholar]
 Lambrechts, M., & Lega, E. 2017, A&A, 606, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., Lega, E., Nelson, R. P., Crida, A., & Morbidelli, A. 2019, A&A, 630, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lega, E., Crida, A., Bitsch, B., & Morbidelli, A. 2014, MNRAS, 440, 683 [Google Scholar]
 Lega, E., Morbidelli, A., Bitsch, B., Crida, A., & Szulágyi, J. 2015, MNRAS, 452, 1717 [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Machida, M. N., Kokubo, E., Inutsuka, S.i., & Matsumoto, T. 2008, ApJ, 685, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Machida, M. N., Kokubo, E., Inutsuka, S.I., & Matsumoto, T. 2010, MNRAS, 405, 1227 [Google Scholar]
 Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marleau, G.D., Kuiper, R., Béthune, W., & Mordasini, C. 2023, ApJ, 952, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (Massachusetts: Courier Corporation) [Google Scholar]
 Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icarus, 232, 266 [Google Scholar]
 Morbidelli, A., Batygin, K., & Lega, E. 2023, A&A, 675, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2 [Google Scholar]
 Nelson, R. P., Lega, E., & Morbidelli, A. 2023, A&A, 670, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [Google Scholar]
 Ronnet, T., & Johansen, A. 2020, A&A, 633, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schulik, M., Johansen, A., Bitsch, B., & Lega, E. 2019, A&A, 632, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schulik, M., Johansen, A., Bitsch, B., Lega, E., & Lambrechts, M. 2020, A&A, 642, A187 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stolker, T., Marleau, G. D., Cugno, G., et al. 2020, A&A, 644, A13 [EDP Sciences] [Google Scholar]
 Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65 [CrossRef] [Google Scholar]
 Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853 [Google Scholar]
 Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Transport and Accretion in PlanetForming Disks (Tucson: University of Arizona Press) [Google Scholar]
 Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, J. J., Ginzburg, S., Ren, B., et al. 2020, AJ, 159, 263 [Google Scholar]
 Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148 [Google Scholar]
 Zhou, Y., Bowler, B. P., Wagner, K. R., et al. 2021, AJ, 161, 244 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Jiang, Y.F., Baehr, H., et al. 2021, MNRAS, 508, 453 [NASA ADS] [CrossRef] [Google Scholar]
except for the small inner bump in the Lm disk where viscous heating dominates; see Fig, A.1, bottom right
All Tables
Total mass inflow, $\left{\dot{M}}_{i}\right$ and outflow, ${\dot{M}}_{o}$ rates.
All Figures
Fig. 1 Radial density profiles (top), azimuthally and vertically averaged temperature (middle), and disk aspect ratio (bottom) profiles for the three models at the end of phase A (left) and C (right). The grey dashed lines indicate the boundary of the Hill sphere and the black dotted line the position of the planet. The density of the Lm and LmLv cases has been multiplied by 10 in order to superpose the density background and appreciate the differences in the gap’s depth and width. We notice that in the LmLv case a large amount of gas is accumulated in a radial domain corresponding to the Hill sphere. This is not gas accumulated around the planet but it is the signature of the libration region around the L4 and L5 Lagrangian points not being depleted yet (a known fact in the case of lowviscosity disks). In the middle panel, the midplane temperature is overplotted with same colours (dashed) on a slice at the planet’s azimuth φ = 0. The higher resolution associated with the deepening of the potential (decreasing smoothing) results in a huge temperature increase up to [1600, 1100, 600] K in the nominal model, the Lm case, and the LmLv model, respectively. In the bottom panel, we remark that it is only in the nominal case that the formation of the gap modifies the aspect ratio considerably. 

In the text 
Fig. 2 Vertically integrated density of the simulated nominal disk (simulation N in Table 1) in the highresolution phase (phase C, see text) with nonuniform grid geometry. We clearly see the spiral arms originating at the planet position (x_{p}, y_{p}) = (5.2, 0) [au], revealing the interaction between the planet and the disk. The disk annulus fully contains the gap and slightly extends over the edges. 

In the text 
Fig. 3 Parameter details for all three models. Left: midplane density field interpolated in cartesian planetocentric coordinates(x, y) = (R cos(ϕ), R sin(ϕ)) in units of the Hill radius). The cartesian components (υ_{x}, υ_{y}) of the midplane velocity field are plotted with arrows. The coloured lines are straightforward (full) and backward (dashed) integration of streamlines obtained with the 2D midplane velocity field. Right: azimuthally averaged density field interpolated in planetocentric coordinates and represented in $(\stackrel{~}{R},z)=(R\mathrm{sin}(\vartheta ),R\mathrm{cos}(\vartheta ))$ plane. The azimuthally averaged and density weighted velocity is represented with arrows. The colours of the arrows correspond to values of the of the gas speed normalised over the sound speed(respectively the midplane speed (left column) and the absolute value of υ_{z} (right column)). The yellow contour levels correspond to the vertical component of the angular momentum in units of the Keplerian momentum with respect to the planet. The coloured markers give the initial conditions for the streamlines of Table 2, projected on the (x, y) and ($\stackrel{~}{R}$, z) planes, respectively, and the black solid lines identify the 0 energy contour (see Eq. (4)). 

In the text 
Fig. 4 Parameter details for all three models. Topleft: radiative diffusion timescale (from Eq. (6)). Topright: opacity. Bottomleft: volume density. Bottom right: temperature. All the quantities are obtained through averages over spherical shells centred on the planet. 

In the text 
Fig. 5 3D representation of streamlines with initial conditions on the Hill sphere (see Table 2). The colour bar corresponds to planetocentric mass flow rate on the Hill surface. The sign is negative for mass inflow and positive for mass outflow. The total integration time is short (about four Jupiter orbits) in order to make this plot clear. The evolution of the same streamlines on a longer integration allows us to follow the dynamics on a projected plane (x/R_{H}, y/R_{H}) and to appreciate the evolution of the distance with respect to the planet in Fig. 6. 

In the text 
Fig. 6 Projection of the streamlines of Fig. 5 on the planetocentric x − y plane (left). The blue contours give the boundaries of the disk. Gas approaching the Hill sphere can have a rapid transit phase in the horseshoe region and be released in the inner or the outer disk trapped in the planetary system. Vertical evolution of streamline with respect to time (middle). Evolution of the distance with respect to the planet as a function of time (right). The dashed lines correspond to backward integration: in the left panel we see that gas approaching the Hill sphere comes from Uturn motion passing close to the planet from the inner or outer branch of the horseshoe region. 

In the text 
Fig. 7 Same as Fig. 3, right panels showing a zoom close to the planetary system. Vectors of components $({v}_{\stackrel{~}{R}},{v}_{z})$ are normalised and the colours correspond to the norm values. The norms cover a range of about five orders of magnitude, with smaller values deep inside the planet’s envelope. 

In the text 
Fig. 8 Absolute value of the inflow and outflow mass rate as a function of z on the Hill sphere. In the Lm and LmLv cases the mass inflow is larger than the mass outflow except very close to the midplane which appears decreating. Instead, in the N case, the Hill sphere is accreting up to z = 0.2[R_{H}] and for z > 0.7[R_{H}]. In all the models, the contribution to the total mass flow in the Hill sphere becomes negligible at z > 0.8[R_{H}]. 

In the text 
Fig. 9 Accretion rate as a function of the distance to the planet in the Hill sphere computed interpolating the time series of the mass contained in a sphere of radius R over 9 planet’s orbital periods (see text). For each case, we report the rate of inflow (circles) and the accretion rate on the Hill sphere (triangle), as computed on a snapshot at five orbital periods (see Table 3). 

In the text 
Fig. A.1 Midplane temperature (top) and disks’ aspect ratio (bottom) as a function of radial distance from the host star for 2D simulations for the nominal and Lm cases (left and right, respectively). The red lines include the effect of stellar irradiation. The other simulations have no stellar heating and differ in the boundary temperature, T_{b}, at disk’s surface. The value of T_{b} that allows us to recover the profile of the stellar irradiated disk is used for 3D simulations of in the rest of the paper. Horizontal dotted lines in the top plots correspond to the values of T_{b} in the simulations. The vertical black dotted line marks the position of a Jupitermass planet in the paper. 

In the text 
Fig. B.1 Same as Fig. 4 but for the Lm case (green curves) and for higher resolution simulations with: a smoothing length ϵ = 0.2 (red curves) and ϵ = 0.1 (blue curves). We remark that decreasing the smoothing length impacts the innermost regions (x < 0.1R_{H}) leaving unchanged the outer part of the Hill sphere. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.