The 3D structure of disc-instability protoplanets

Context. The model of disc fragmentation due to gravitational instabilities offers an alternate formation mechanism for gas giant planets, especially those on wide orbits. Aims. Our goal is to determine the 3D structure of disc-instability protoplanets and to examine how this relates to the thermal physics of the fragmentation process. Methods. We modelled the fragmentation of gravitationally unstable discs using the SPH code PHANTOM, and followed the evolution of the protoplanets formed through the first and second-hydrostatic core phases (up to densities 1e-3 g/cm3). Results. We find that the 3D structure of disc-instability protoplanets is affected by the disc environment and the formation history of each protoplanet (e.g. interactions with spiral arms, mergers). The large majority of the protoplanets that form in the simulations are oblate spheroids rather than spherical, and they accrete faster from their poles. Conclusions. The 3D structure of disc-instability protoplanets is expected to affect their observed properties and should be taken into account when interpreting observations of protoplanets embedded in their parent discs.

Gravitational instabilities develop in protostellar discs when the Toomre criterion (Toomre 1964) is satisfied, where c s is the sound speed, κ is the epicyclic frequency, and Σ is the surface density of the disc, at an orbital radius R. Gravitational instability leads to disc fragmentation when the disc cools sufficiently quickly, t cool < (0.5−2)t orb (i.e. a few orbital periods).Magnetic fields may also play an important role in disc formation (Wurster & Li 2018;Lebreuilly et al. 2024;Hennebelle et al. 2020) and subsequent disc fragmentation (Commerçon et al. 2010).It is believed that magnetic fields tend to act towards suppressing disc fragmentation, although this may still happen under the appropriate conditions (Commerçon et al. 2010;Forgan et al. 2017;Deng et al. 2021).
The fragments produced by gravitational instability have masses that are a few times the mass of Jupiter (M J ), but the final mass they acquire may be much higher (Stamatellos & Whitworth 2009a;Kratter et al. 2010;Vorobyov 2013;Kratter & Lodato 2016;Mercer & Stamatellos 2017;Fletcher et al. 2019).The disc instability theory naturally forms gas giant planets on wide orbits (Stamatellos & Whitworth 2009a), where both criteria for disc fragmentation are satisfied.However, interactions with passing stars may destroy an initial population of such planets (Carter & Stamatellos 2023), in line with direct imaging observations (e.g.Bowler & Nielsen 2018;Vigan et al. 2021) that show that massive gas giants on wide orbits are not very common (only a small percentage of stars host such planets, up to a maximum of 5-10% of stars, with a small dependence on the stellar host mass).
The evolution of disc-instability fragments to protoplanets goes though the phases of the first and second hydrostatic cores (Larson 1969;Masunaga & Inutsuka 2000;Stamatellos et al. 2007b;Stamatellos & Whitworth 2009b).This means that the initial stages of the evolution of disc-instability planets are similar to those of a star within a collapsing molecular cloud core, albeit at a much smaller scale: the initial core (i.e. the fragment formed by disc instability that will evolve to a planet) has a size of a few AU, has a mass of a few M J , and is rather rapidly rotating (Stamatellos & Whitworth 2009b;Mercer & Stamatellos 2020).
Observations of planets still embedded in their parent discs (usually referred to as protoplanets) have become possible in the last few years (see review by Currie et al. 2023).The two protoplanets around the 5 Myr old star PDS 70 are the first unambiguous discoveries (Keppler et al. 2018;Haffert et al. 2019).They orbit at distances of 20 and 34 AU from the central star.PDS 70 b has an estimated mass of <12 M J ; the mass of PDS 70 c is uncertain.These protoplanets show signatures of gas accretion, as evidenced by H α emission (Wagner et al. 2018;Haffert et al. 2019), and are attended by circumplanetary discs (Stolker et al. 2020;Benisty et al. 2021).Recently, a protoplanet has been discovered around Aurigae AB (Currie et al. 2022), a 1−3 Myr old star.This protoplanet has an estimated mass of ∼9 M J and orbits at ∼93 AU from its parent star.
As more direct (and indirect) observations of protoplanets are likely in the near future, it is important to determine their properties when they form by different scenarios (core accretion and disc instability) so that we may identify the dominant gas giant planet formation mechanism.In this work we perform a set of hydrodynamic simulations of disc fragmentation to determine the 3D structure of disc-instability protoplanets.In Sect. 2 we discuss the disc initial conditions and the methods used for the simulations.In Sect. 3 we present their general results, and in Sect. 4 the density, temperature, and velocity profiles of the protoplanets that form in the simulations.In Sect. 5 we focus on the shape of disc-instability protoplanets, and in Sect.6 we summarize the main results of our study.

Methodology
We model the thermodynamics of gravitationally unstable discs with the smoothed particle hydrodynamics code PHANTOM (Price et al. 2018), using a barotropic equation of state (e.g.Bate 1998).We vary the density at which the equation of state switches from isothermal to adiabatic, the adiabatic index, and the initial disc temperature as this is set by stellar heating.

Disc initial conditions
We set a disc with mass of M D = 0.6 M around a host star of 0.8 M .The disc extends from 10−300 AU and it is represented by N SPH = 4 × 10 6 particles.The disc mass is chosen so that many fragments can form due to disc fragmentation in each simulation and facilitate a statistical study of their properties.The minimum mass that can be resolved is N neigh M D /N SPH ∼ 7.5 × 10 −4 M J , which is much lower than the opacity limit for fragmentation (a few M J ; e.g.Whitworth & Stamatellos 2006).Therefore, disc fragmentation is property resolved.
The surface density of the disc is set to where R in = 10 AU is the inner disc radius, and Σ 0 = 1.53 × 10 3 g cm −2 .The disc temperature profile is set to where T 1 AU = [150, 200] K.The above disc initial conditions ensure that the disc is initially Toomre unstable outside ∼50 AU.

Disc thermodynamics
Hydrodynamic simulations often use a barotropic equation of state (i.e.P ∝ ρ γ ) to reproduce the results of more computationally exhaustive radiative hydrodynamic simulations (Masunaga et al. 1998;Masunaga & Inutsuka 2000;Whitehouse & Bate 2004;Mercer et al. 2018).To emulate the thermal effects during gravitational fragmentation in protostellar discs we use a hybrid four-piece barotropic equation of state that is modified to include radiative feedback from the central protostar.More specifically the temperature of an SPH particle i is ).The disc becomes gravitationally unstable and fragments.Four of the fragments or protoplanets are followed until they reach density 10 −3 g cm −3 .
where T (R i ) is set by the central star (see Eq. ( 3)) and T B (ρ i ) is provided from the barotropic equation, where γ 1 , γ 2 , and γ 3 are the adiabatic indices that control the stiffness of the equation of state in the three density regions (i.e.how fast the gas temperature rises due to compressional heating during the collapse).The first region (ρ < ρ 1 ; typically ρ 1 ∼ 10 −13 g cm −3 , T 10 K) corresponds to the phase of isothermal collapse where the gas is optically thin and its radiation escapes freely.The second region (ρ 1 < ρ < ρ 2 ; typically ρ 2 ∼ 3 × 10 −12 g cm −3 , T ∼ 10−100 K, γ = 5/3) corresponds to the phase where the gas is optically thick and starts heating up.The third region (ρ 2 < ρ < ρ 3 ; typically ρ 3 ∼ 6 × 10 −9 g cm −3 , T ∼ 100−2000 K, γ = 7/5) corresponds to the phase where the rotational degrees of molecular hydrogen have been excited.Finally, the last region (ρ > ρ 3 ; T > 2000 K, γ = 1.1) corresponds to the phase where molecular hydrogen starts to dissociate.
The first critical density, ρ 1 , effectively determines when the fragment becomes optically thick, and is therefore a measure of of the disc opacity and metallicity.The critical densities ρ 2 and ρ 3 are set (for each γ 1 , γ 2 ) so as to correspond to temperatures of 100 K and 2000 K, where the rotational degrees of molecular hydrogen are excited and the dissociation of molecular hydrogen commences, respectively.

Disc fragmentation and protoplanet formation
We performed high-resolution disc fragmentation simulations with nine different sets of parameters.The parameter sets investigated are summarized in Table A.1.The disc initial conditions were chosen so that the discs quickly become gravitationally unstable, as evidenced by strong spiral arms, and fragment.These self-gravitating fragments are referred to as protoplanets.We followed their evolution to density 10 −3 g cm −3 .Simulations L6, page 2 of 9 with stiffer equations of state (γ = 1.66) form fewer protoplanets, due increased compressional heating providing support against collapse.A total of 107 protoplanets form in all simulations.
A typical outcome of a simulation (the benchmark simulation) is shown in Fig. 1.Four of the protoplanets that form in the simulation are followed until they reach a density of 10 −3 g cm −3 .Two of these protoplanets are shown in more detail in Fig. 2, where they are plotted at two different times.We also plot representative protoplanets, as seen face-on and edge-on, for each set of parameters (see Figs. A.1-A.4).The non-axisymmetric and flattened morphology of these protoplanets is evident.
The evolution of a fragment or protoplanet once it starts collapsing follows the same stages as the collapse of a solarmass molecular cloud to a protostar (Stamatellos & Whitworth 2009b).The main difference is that the mass of the fragment itself is much lower, of the order of 10−100 M J .The collapse of the fragment or protoplanet is initially isothermal, with the temperature set by how far away the protoplanet is from the central star (typically 10−30 K).Once the protoplanet becomes optically thick, the first hydrostatic core forms (Larson 1969;Stamatellos et al. 2007b), which grows in mass, and slowly contracts and heats up; an accretion shock forms around the first core as infalling gas decelerates.When the temperature rises to 2000 K the molecular hydrogen dissociation initiates the second collapse and the second hydrostatic core forms (see Stamatellos & Whitworth 2009b;Mercer & Stamatellos 2020).The mass of the first core is of the order of 10−20 M J , whereas the mass of the second core is a few M J .The final mass of the protoplanet will be decided through interactions with the disc (Mercer & Stamatellos 2020).Each protoplanet is represented by at least 6 × 10 5 SPH particles, and therefore the thermodynamics of the collapse is properly resolved (Stamatellos et al. 2007b).The property of the protoplanets discussed later on refers to those when the density of 10 −3 g cm −3 is reached at their centres.

Three-dimensional structure of disc-instability protoplanets
Previous studies (e.g.Mercer & Stamatellos 2020) have assumed that disc-instability protoplanets, are spherically symmetric.However, these protoplanets form in a disc with a nearly Keplerian rotational profile, such that they are rotating.Therefore, they are expected to be flattened, as happens in rotating collapsing clouds leading to protostar formation (Bate 1998;Saigo & Tomisaka 2006;Saigo et al. 2008).
To investigate the 3D structure of the protoplanets formed in the simulations in more detail, we calculate the density, temperature, rotational velocity, and infall velocity along different directions from the centre (±x, ±y, ±z) of each protoplanet .We also calculate axisymmetric averages on the x−y plane (which is assumed to be the plane of rotation) and spherical averages.The results for a typical protoplanet are shown in Fig. 3 (for the protoplanet formed in the benchmark simulation, see Fig. 2a).The density in the z-direction drops faster with radius than the density in the other two directions indicating that the protoplanet is flattened.The densities in the ±x and ±y directions are very similar apart from the edges of the protoplanet.This is also true for the temperature profile of the protoplanet.The rotational velocity of the protoplanet (which is not calculated in the z-direction, as this is the axis of rotation) shows differences along different directions as a result of the formation environment of the protoplanet that is being fed with gas from the disc.The infall velocity (see Fig. 3) is considerably higher along the poles of the protoplanet (i.e. in the ±z directions).The presence of the accretion shocks L6, page 3 of 9 around the first and second core is evidenced by the maxima and minima in the infall velocities (see Fig. 3), which correspond to the start of the shock when gas starts to decelerate before it falls onto the first and second core, respectively, and stops.
In Fig. 4 we compare the structure of the four protoplanets that form in the benchmark run.We plot axisymmetric averages (solid lines), and averages along the ±z direction (dotted lines) for the density, temperature, rotational velocity, and infall velocity.The four protoplanets show similar density and temperature profiles apart from the edges, where they interact with the protostellar disc.However, the rotational velocity and infall velocity of gas towards their centres show significant differences that are indicative of their different formation histories.In all cases the infall velocity along the poles (i.e. in the ±z directions) of the protoplanet is much higher (a factor of ∼2) than that along the protoplanet equator (i.e. on the x−y plane).This is demonstrated in Fig. 5 where the velocity vectors of the gas are plotted over the density at different planes for the protoplanet in Fig. 2, top.We find that accretion of gas onto protoplanets happens from the polar directions, as is also expected for core-accretion planets (Tanigawa et al. 2012).
Our results for all protoplanets that form in the simulations show that (i) protoplanets are flattened and symmetric with respect to the x−y plane (i.e. the protostellar disc midplane), and (ii) protoplanets are nearly axisymmetric, although there are differences near their edges due their formation environment (e.g.interactions with spiral arms) and formation history (i.e. when a protoplanet forms due to a collision between two fragments) (see ).For simplicity, we assume for the rest of the discussion that protoplanets are axisymmetric so that we can make comparisons between protoplanets formed in different simulations.We note that our study has neglected the effect of magnetic fields, which may influence disc formation and subsequent fragmentation (e.g.Commerçon et al. 2010), and therefore the 3D structure of protoplanets.

The shape of disc-instability protoplanets
Disc-instability protoplanets are nearly axisymmetric (with respect to the rotation axis z) and symmetric with respect to the x−y plane; therefore, they can be described as oblate spheroids.).Gas infall velocities onto the protoplanet are asymmetric; higher velocities are seen towards the poles of the protoplanet.
To quantify their shape we use three metrics.The first metric is the aspect ratio of the first core, e fc , as this is calculated at the inner boundary of the accretion shock around it (i.e.where the infall velocity is minimum or almost zero).This is the ratio of the inner first core radius on the x−y plane over the corresponding radius in the z direction (see Fig. 3).This ratio cannot be calculated accurately for all protoplanets as the outer region of the protoplanet is of low-density and therefore not L6, page 4 of 9 .Aspect ratios of the first core e fc , the second core, e sc , and the fiducial core, e ρ , (top to bottom) of all protoplanets that form in the simulations.The colour of each symbol corresponds to different adiabatic indices, the shape to different density threshold for the gas becoming optically thick (i.e. to different disc opacities and metallicity), and the type of each symbol (filled or unfilled) to different stellar radiation fields (as marked on the graph legend).
well-resolved in some cases.The second is the aspect ratio of the second core, e sc , as this is calculated at the outer boundary of the accretion shock around it (i.e.where the infall velocity is maximum).The third metric is the aspect ratio, e ρ , using a fiducial distance around the centre of the protoplanet where the density drops to ρ c = 10 −9 g cm −3 (which roughly corresponds to the radius of the first core), i.e. the ratio of the fiducial radius in the z direction over the fiducial radius on the x−y plane.This metric has the advantage that can be defined for all protoplanets.
All three metrics paint the same picture regarding the protoplanets' morphology (see Fig 6).Most protoplanets have aspect ratios <1, which means they are flattened, oblate spheroids rather than spherically symmetric.A few protoplanets have very high aspect ratios (above 1); these are possible outcomes of merger events.
In Fig. 7 we plot the aspect ratios of the first and second core with respect to the corresponding ratios of the rotational to gravitational energy (β fc and β sc , respectively).We see that second cores with higher β sc values tend to be flat, but there is no such relation for the first cores.This suggests that the shape of the first cores is determined by interaction with the disc, whereas the shape of the second cores is due to their rotation.However, there are a few cases with high second core aspect ratio even with high rotational-to-gravitational energy ratios, suggesting a violent formation process (e.g.mergers or strong interactions with spiral arms).A comparison of the distributions of the aspect ratios of the first and second cores is shown in Fig. 8.The two distributions are distinctly different, with second cores aspect ratios peaking around e sc ∼ 0.7−1, whereas first cores are highly flattened, with aspect ratios peaking around e fc ∼ 0.1, which is similar to the disc scale height.
A stiffer equation of state (γ 1 = 1.66 vs. γ 1 = 1.4) generally results in more spherical first cores ( e fc = 0.62 vs. e fc = 0.26), but flatter second cores ( e sc = 0.68 vs. e sc = 0.96; compare the different colours in Fig. 6).This occurs because a stiffer equation of state S also results in second cores with higher rotational-to-gravitational energy ratios ( β sc = 0.27 vs. β sc = 0.17; see Fig. 7).For the first cores the β-ratios are similar β fc = 0.23.
A disc with higher metallicity and opacity (ρ 1 = 6 × 10 −13 g cm −3 vs. ρ 1 = 10 −13 g cm −3 ) results in slower rotating ( β fc = 0.19 vs. β fc = 0.28), but slightly more spherical first cores ( e fc = 0.41 vs. e fc = 0.38), and in slower rotating ( β sc = 0.17 vs. β sc = 0.21, but slightly flatter second cores ( e sc = 0.84 vs. e sc = 0.99).However, these differences are rather small (compare triangles with circles in Figs. 6 and 7).A disc with greater heating from the central star (T 1 AU = 200 K vs. T 1 AU = 150 K) has no effect on the rotational-to-gravitational L6, page 5 of 9 energy ratios, but the first and second cores are slightly more spherical ( e fc = 0.46 vs. e fc = 0.33, e sc = 0.92 vs. e sc = 0.87, respectively) (compare filled with unfilled symbols in Figs. 6 and 7)).Finally, there seems to be no correlation between the shapes of the first and second cores and the position where they form within the disc.

Conclusions
Disc-instability protoplanets are not spherically symmetric, but close to being oblate spheroids.Their outer regions show more complex, asymmetric structure due to interactions with the protostellar disc and their formation history.Gas accretion happens faster from the protoplanet poles than from the protoplanet equator.We expect that this may lead to a strong modification of the observed properties of protoplanets (e.g.their spectrum, H α emission; see Zhu 2015;Marleau et al. 2022Marleau et al. , 2023) ) with the viewing angle that needs to be taken into account when interpreting observations, like those of PDS 70 b,c (Keppler et al. 2018;Haffert et al. 2019) and AB Aurigae b (Currie et al. 2022).

Appendix A: Simulation parameters and gallery of protoplanet surface density plots
The parameter sets investigated are listed in Table A.1.For an explanation of the different parameters see Section 2.
Representative protoplanets for each of the eight sets of equation of state parameters are plotted in Figs.A.1-A.4.The different morphologies of the protoplanets for different equation of state parameters are clearly seen in these plots.

Fig. 1 .
Fig. 1.Surface density of the benchmark run disc (in g cm −2).The disc becomes gravitationally unstable and fragments.Four of the fragments or protoplanets are followed until they reach density 10 −3 g cm −3 .

Fig. 2 .
Fig. 2. Surface density of two of the protoplanets (one per row; face-on view) that form in the benchmark run (in g cm −2 ), plotted when their central density is ρ c = 10 −9 (left column) and ρ c = 10 −5 g cm −3 (right column).

Fig. 3 .
Fig. 3. Density, temperature, rotational velocity, and infall velocity at different directions from the centre of one of the protoplanets that form in the benchmark run (see Fig. 2, top) as marked on the graph.The axisymmetric averages are represented by the black dotted lines and the spherical averages are shown by the black dashed line.

Fig. 4 .
Fig. 4. Radial profiles on the x − y plane (solid lines) and vertical profiles (±z direction; dotted lines) of the density, temperature, rotational, and infall velocity (panels a, b, c, and d, respectively) for the four protoplanets that form in the benchmark simulation.

Fig. 5 .
Fig. 5. Velocity vectors of the gas flow onto a disc-instability protoplanet (see Fig. 2, top) on the x−y plane (top), x−z plane (middle), and y−z plane (bottom), overplotted on the corresponding densities (in units of g cm −3).Gas infall velocities onto the protoplanet are asymmetric; higher velocities are seen towards the poles of the protoplanet.
Fig.6.Aspect ratios of the first core e fc , the second core, e sc , and the fiducial core, e ρ , (top to bottom) of all protoplanets that form in the simulations.The colour of each symbol corresponds to different adiabatic indices, the shape to different density threshold for the gas becoming optically thick (i.e. to different disc opacities and metallicity), and the type of each symbol (filled or unfilled) to different stellar radiation fields (as marked on the graph legend).

Fig. 7 .Fig. 8 .
Fig. 7. Aspect ratios of the first and second core with respect to the ratios β fc , and β sc , of the rotational to the gravitational energy of the first and second core, respectively.Symbols as in Fig.6.Second cores are generally more flattened when they rotate faster, but first cores do not show such dependence.

Fig. A. 3 .
Fig. A.3.Same as in Fig. A.1, but for protoplanets in Runs 5 − 8. Fig. A.4. Same as in Fig. A.2, but for protoplanets in Runs 5 − 8. Due to the projection, more than one protoplanet appears in some of these plots.

Table A .
1. Equation of state parameters used for the disc fragmentation simulations (see Section 2).