Identifying and Analysing Protostellar Disc Fragments in Smoothed Particle Hydrodynamics Simulations

TL;DR

This paper introduces a new density derivative method for identifying protostellar disc fragments in SPH simulations, compares it with an existing method, and analyzes fragment survival and dynamics.

Contribution

The study develops a novel density derivative technique for detecting disc fragments and demonstrates its effectiveness alongside the CLUMPFIND method in SPH simulations.

Findings

All fragments are detected using the density derivative method.

Fragments that survive are only detected by the gravitational potential method.

A tentative relationship between azimuthal wavenumber and maximum semi-major axis is proposed.

Abstract

We present a new method of identifying protostellar disc fragments in a simulation based on density derivatives, and analyse our data using this and the existing CLUMPFIND method, which is based on an ordered search over all particles in gravitational potential energy. Using smoothed particle hydrodynamics, we carry out 9 simulations of a $0.25$ M$_{\odot}$ disc around a 1 M$_{\odot}$ star, all of which fragment to form at least 2 bound objects. We find that when using all particles ordered in gravitational potential space, only fragments that survive the duration of the simulation are detected. When we use the density derivative method, all fragments are detected, so the two methods are complementary, as using the two methods together allows us to identify all fragments, and to then determine those that are likely to be destroyed. We find a tentative empirical relationship between the dominant azimuthal wavenumber in the disc $m$ and the maximum semi-major axis a fragment may achieve in a simulation, such that $a_{\mathrm{max}}\propto\frac{1}{m}$. We find the fragment destruction rate to be around half that predicted from population synthesis models. This is due to fragment-fragment interactions in the early gas phase of the disc, which can cause scattering and eccentricity pumping on short timescales, and affects the fragment's internal structure. We therefore caution that measurements of eccentricity as a function of semi-major axis may not necessarily constrain the formation mechanism of giant planets and brown dwarfs.

Tables

Table 1. Table 1 : Parameter values of the Gaussian fits applied to the histograms in Figure 13 . From left to right, the columns are plot, mean, standard deviation, amplitude (without offset) and offset constant.

Plot	$\mu$	$\sigma$	$A$	$k$
Initial ecc.	0.094	0.095	0.070	0.006
Final ecc.	0.107	0.071	0.079	0.006
Initial inc. (no outlier)	0.052	0.029	0.173	0.003
Initial inc. (with outlier)	0.053	0.028	0.051	0.006
Final inc. (no outlier)	0.001	0.001	0.119	0.005
Final inc. (with outlier)	0.001	0.001	0.063	0.007

Equations

$$\mathcal{P}=0.03\times 10^{2.0[\mathrm{Fe}/\mathrm{H}]},$$

$$\frac{\partial {\rho }_i}{\partial r}=\frac{1}{n}\sum _{j=-\frac{n-1}{2}}^{j=\frac{n-1}{2}}\frac{\partial {\rho }_{i+j}}{\partial r}$$

$$\frac{\partial {\rho }_i}{\partial \varphi }=\frac{1}{n}\sum _{j=-\frac{n-1}{2}}^{j=\frac{n-1}{2}}\frac{\partial {\rho }_{i+j}}{\partial \varphi }$$

$$A_{\mathrm{m}}=\Bigg{|}\sum_{i=1}^{N_{\mathrm{region}}}\frac{e^{-im\phi_{i}}}{N_{\mathrm{region}}}\Bigg{|},$$

$$a_{\max }\propto \frac{1}{m}.$$

$$i=\arccos\Bigg{(}\frac{L_{\mathrm{z}}}{|\vec{L}|}\Bigg{)},$$

$$F_{\mathrm{p}}=A{e}^{-\frac{1}{2}(\frac{x-\mu }{\sigma }{)}^2}+k,$$

$$P=K{\rho }^{\frac{n+1}{n}},$$

$$\frac{|B|}{t_{\mathrm{ff}}}=\Bigg{(}\frac{3}{5\pi}\frac{\mathrm{G}^{3}M^{5}}{R^{5}}\Bigg{)}^{\frac{1}{2}},$$

$$B=-\frac{3}{5}\frac{\mathrm{G}{M}^2}{R},$$

$$t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32\mathrm{G}\rho }}.$$

$$\frac{|B|}{t_{\mathrm{ff}}}\lesssim L,$$

$$L=e4\pi R^2\sigma T^4,$$

$$M_{crit}=\Bigg{(}\frac{400}{9}\frac{e^{2}\pi^{4}\sigma^{4}T^{8}R^{9}}{\mathrm{G}}\Bigg{)}^{\frac{1}{5}}$$

$$M_{\mathrm{J}}\lesssim M\lesssim M_{\mathrm{crit}},$$

$$M_{\mathrm{frag}}\approx 6\frac{T^{\frac{1}{4}}}{e^{\frac{1}{2}}}{\mathrm{M}}_{\mathrm{J}}.$$

$${\Gamma }_{\mathrm{G}}=\mathrm{sgn}(k)\frac{{\pi }^{2}m{R}_0\mathrm{G}{\Sigma }_1^2}{k^2},$$

$$k\equiv \frac{\partial (mf(R))}{\partial R},$$

$${\Gamma }_{\mathrm{G}}\propto \frac{1}{m}.$$

$${\Gamma }_{\mathrm{G}}$$

$$\displaystyle=\frac{d}{dt}\bigg{(}m|\vec{v}\times\vec{r}|\bigg{)}$$

$$\displaystyle=\frac{d}{dt}\bigg{(}m|\vec{v}||\vec{r}|\mathrm{sin}[\theta]\bigg{)},$$

Full text

††pagerange: Identifying and Analysing Protostellar Disc Fragments in Smoothed Particle Hydrodynamics Simulations–Identifying and Analysing Protostellar Disc Fragments in Smoothed Particle Hydrodynamics Simulations††pubyear: 2016

1 Introduction

There are two distinct modes of planet formation in protostellar discs. The first, and most widely accepted, is the core accretion model (CA) (Pollack et al., 1996; Hubickyj et al., 2005). In this model, growth begins with dust grains of $\sim 1\/\micron$ that coagulate rapidly into larger particles, ultimately settling into the disc midplane where there is enough material for them to grow to kilometre-sized planetesimals. These planetesimals can then grow via collisions into planetary cores, and if sufficiently massive, and if the gas disc has not dissipated, will accrete a gaseous envelope, ultimately becoming a gas giant planet (Pollack et al., 1996; Lissauer, 1993).

Most observational evidence favours this formation mechanism. For example, gas giant planets are preferentially found around metal-rich stars (Santos et al., 2004), with an empirical relationship that quantifies the probability, $\mathcal{P}$, of gas giant planet formation as

[TABLE]

where $[\mathrm{Fe}/\mathrm{H}]$ is the metallicity of the host star relative to solar metallicity (Fischer & Valenti, 2005). Numerical work (Cai et al., 2005) has suggested that this would not be the case if the second mode of planet formation, gravitational instability (GI), were the dominant formation mechanism of these planets, since an increase in metallicity responds to a decrease in cooling rate, resulting in weaker GI activity. This ultimately decreases the likelihood of these systems fragmenting, since weak GI corresponds to smaller stresses in the disc. On the other hand, it has also been shown that metallicity variation makes very little difference to the occurrence of fragmentation (Boss, 2002).

In the GI scenario, gas giant planets and brown dwarfs form by direct gravitational collapse in the gaseous protostellar disc (Kuiper, 1951; Cameron, 1978; Boss, 1997, 1998). This happens rapidly, in a relatively early phase of the disc’s life when it is massive enough to be self-gravitating. The advantage of this mechanism is its rapidity; gas giants are able to form on timescales shorter than typical disc dispersion timescales ($\sim$ 5 Myr; Haisch et al. 2001). While CA is certainly the most widely accepted model, there are barriers to grain growth at several length scales which seem to indicate difficulty in forming planetary mass objects within the disc lifetime. The most famous of these is the so-called metre barrier; as grains increase in size, so do their relative velocities, which makes grain fragmentation, rather than coagulation, the most likely outcome.

A promising solution to this problem is the pebble accretion theory (Lambrechts & Johansen, 2012; Levison et al., 2015). Pebbles are grouped together due to the streaming instability (Youdin & Goodman, 2005), whereby solid particles orbit at Keplerian velocity, but the gas is pressure supported from the host stellar radiation, causing the gas to orbit at sub-Keplerian speeds. Feeling a headwind, solids slow, losing angular momentum and migrating inwards. As more solid particles migrate inwards, they cluster together, and if the solid-to-gas ratio is sufficiently large (order unity Youdin & Goodman 2005), then the backreaction of the dust on the gas will change the local gas velocity. This, in turn, alters the local drag force, promoting the pile up of solids which may gravitationally collapse if they become sufficiently large (Youdin, 2011). These groups of solid particles may begin to accrete pebbles until they form giant planet cores (Lambrechts & Johansen, 2012).

At smaller scales, the bouncing barrier prevents coagulation of dust grains, as particles of a given size, above a certain velocity, are more likely to bounce off each other than they are to coagulate. This results in growth typically halting at around $\sim$ 1 mm in size. However, there is evidence to suggest that this could be beneficial to planetesimal formation, since the introduction of a few $\sim$ cm sized grains (e.g, through radial drift) can act as catalyst to grain growth, sweeping up grains, while preventing the growth of too many larger objects which would otherwise smash each other apart (Windmark et al., 2012).

It is generally accepted that disc fragmentation is very unlikely in the inner regions (< 50 au) of a protostellar disc (Rafikov, 2005). However, the outer regions of protostellar discs may well be susceptible to fragmentation, offering a formation mechanism for directly imaged planets such as those in the HR8799 system (Marois et al., 2008; Nero & Bjorkman, 2009; Kratter et al., 2010). Core accretion models struggle to explain objects such as those in HR8799, with four planets orbiting at 14, 24, 34 and 68 au, with masses of $\sim$ 5 MJ (Marois et al., 2008; Marois et al., 2010), since there is not thought to be enough material to form such massive objects at these distances. Additionally, the growth timescales of such objects, through core accretion, typically exceed disc lifetimes by a factor of at least $\sim$ 3, using conservative estimates (Pollack et al., 1996). Gravitational instability may, perhaps, offer an explanation as to the formation mechanism of these systems.

However, it has been suggested that disc fragmentation rarely forms planetary mass objects (Rice et al., 2015), with some hydrodynamics simulations (Stamatellos & Whitworth, 2009) suggesting objects formed by this mechanism quickly grow to brown dwarf masses ($M>$ 13 MJup), with lower limits placed on the fragment mass of $\sim 3-5$ MJ (Kratter et al., 2010; Forgan & Rice, 2011). This is compounded by the recent possible observation, for the first time, of disc fragmentation in action (Tobin et al., 2016), which shows the birth of three protostars that are well above the upper limit of the planetary mass regime.

The recent reformulation of the GI scenario in to what is now known as "tidal downsizing" (Boley et al. 2010; Boley et al. 2011; Nayakshin 2010, 2011a, 2011b) does, however, have positive implications for producing low-mass planets at low semi-major axes. The key is to consider the subsequent evolution of fragments into planetary embryos, through dust growth, radial migration and tidal disruption. Forgan & Rice (2013) combined the physical processes of tidal downsizing with semi-analytic models of disc evolution (Rice & Armitage, 2009) to produce the first population synthesis model for planets formed through GI. Given the similarities between these fragments and “first cores” (see, e.g. Masunaga et al. 1998), they were modelled as polytropic spheres, with polytropic index $n=1.5$.

They found that $\sim 40\%$ of fragments that formed are ultimately tidally destroyed by the central star, and of those that survive $\sim 40\%$ are gas-giant planets with solid cores of $5\sim 10$ earth masses, and the rest are brown dwarfs with no solid core. They also found that low mass embryos tend to remain at larger semi-major axes due to the tidal downsizing process. Out of over 1 million fragments, there was only one terrestrial type planet (core with no gaseous envelope). These results are inconsistent with GI being the dominant planet formation mechanism, but they are certainly consistent with GI forming brown dwarfs and gas giant planets at large radii.

Population synthesis models, by necessity, make simplifying assumptions about the physics that governs the evolution of each planetary system. In particular, interactions between forming protoplanets, and the interaction of the disc with these protoplanets, are not included in the population synthesis models of Forgan & Rice (2013) that we discuss here. In fact, at the time of writing, these effects are not included in any GI population synthesis models. How important these interactions are in determining the final orbital configuration of a system is something that should be carefully considered before further developments are made to such a model. Quantifying the importance of these interactions is difficult, but some headway can be made by performing SPH simulations of fragmenting protostellar discs, and carefully tracking the evolution of fragments’ orbital and physical properties throughout the duration of the simulation.

In this work, we analyse a suite of SPH simulations of fragmenting protostellar discs, identifying fragments using two different methods. The first, based on the CLUMPFIND algorithm (Williams et al., 1994; Smith et al., 2008), is done using the gravitational potential, and the second is a new method that uses density derivatives. We do not use sink particles (Bate et al., 1995) in our simulations, as by using only SPH particles, we are able to determine the fragment internal structure as it migrates through the disc, and better understand the orbital evolution of the fragment, which is sensitive to its radial mass distribution. We discuss the relative merits of our different detection methods for our simulations, which show a variety of fragmentation scenarios. We consider the implications of our results for current GI population synthesis models, and finally we consider the orbital and physical properties of the fragments in the simulations, comparing them to the orbital and physical properties of the population synthesis models of Forgan & Rice (2013).

The paper is ordered as follows: In section 2, we describe our overall method, outlining our chosen formalism of SPH in section 2.1 and detailing the simulation setup in section 2.2. We present our algorithms in section 2.3, describing our new approach in section 2.3.2 and our adaptation of an existing approach in section 2.3.1. We describe our results in section 3, detailing the relative merits of the different approaches in section 3.1. We compare our results to current gravitational instability population synthesis models in section 3.2. We outline orbital and spin properties of our fragments in section 3.3 and 3.4 respectively. We describe density and temperature profiles of fragments with the most interesting histories in section 3.5, and discuss our findings and conclude in section 4.

2 Method

We run a suite of 9 smoothed particle hydrodynamics (SPH) simulations, with the radiative cooling method of Forgan et al. (2009). Each disc uses 4 million particles, fragments to form at least two bound objects, and - aside from the random number seed used to initialise each disc - has identical initial conditions. Since our aim is to take advantage of running a fully hydrodynamic simulation by tracing the evolution of fragment radial profiles and mass distribution, we run the simulations for as long as it is computationally feasible using hydrodynamics only (without switching to sink particles). In practice, this means running each simulation until the most dense clump becomes too dense to calculate the next timestep. We give a brief description of SPH, and our selected radiative transfer formalism, in section 2.1. We discuss our methods of finding fragments in section 2.3.

2.1 Smoothed Particle Hydrodynamics

Smoothed Particle Hydrodynamics (SPH) is a Lagrangian hydrodynamics formalism, that evolves a fluid by means of a distribution of pseudo-particles (Lucy, 1977; Gingold & Monaghan, 1977). There are many review articles about SPH (see e.g. Monaghan 1992, 2005; Rosswog 2009), but the basic idea is that each particle has a position, mass, internal energy and velocity, and these parameters can be interpolated over to give fluid quantities at any position. Density is calculated by interpolation over the mass distribution, and pressure is determined using an equation of state with internal energy. Gravitational forces are usually computed using a TREE algorithm (Barnes & Hut, 1989), and then the discretised energy and momentum equations are solved. Particle velocities are updated using pressure and gravitational forces, and positions are then updated using these velocities. Internal energy changes are computed by calculating $P\mathrm{d}V$ work, viscous dissipation and radiative cooling and heat conduction.

Cooling calculations in SPH are no simple task. Accounting for polychromatic radiative transfer within a hydrodynamics simulation is not possible with current computational resources, and even post-processing a single snapshot with radiative transfer is computationally expensive (Stamatellos & Whitworth, 2005). Historically, approximations to individual features of radiative transfer were used, such as the cooling time formalism: $\dot{u}=-u/t_{\mathrm{cool}}$ (Rice et al., 2003a). Although this parameterisation is useful, allowing us to probe the effects of different cooling timescales in protoplanetary discs, it is somewhat limited, as it only allows us to model energy lost from an SPH particle. Realistically, if energy is lost from one SPH particle, at least some of that energy will be gained by its surrounding neighbours - this is known as radiative transfer.

Since our aim here is to trace the orbital and profile evolution of fragments within protostellar discs, we wish to capture the effects of radiative transfer as far as is feasible. Therefore, the cooling we implement is the hybrid method of Forgan et al. (2009). The details of the algorithm are given in Forgan et al. (2009), however the basic ideas merge the polytropic cooling method of Stamatellos et al. (2007) with the flux-limited diffusion method of Mayer et al. (2007), which builds on conduction modelling work by Cleary & Monaghan (1999) and the flux-limiter described in Bodenheimer et al. (1990). The biggest advantage is the complementary nature of these two methods, energy loss is handled by polytropic cooling (which flux-limited diffusion cannot do), and positive energy exchange between neighbouring particles is handled by flux-limited diffusion (which polytropic cooling cannot do). Since each method handles a separate process, there is no “double counting” in any part of the system’s overall energy, and these separate parts can simply be summed to calculate the total energy change.

2.2 Simulation setup

We run a total of 9 simulations of 0.25 M⊙ discs, with a 1 M⊙ central star, an inner radius of 10 au, an outer radius of 100 au, and a radial density profile of $\Sigma\propto r^{-1}$. All discs are initially identical in global properties, varying only the random number seed (the integer used to set the starting point for a sequence of random numbers) used to initialise the disc. The SPH particles are randomly distributed in $\phi$, where $\phi$ is azimuthal angle, and the $r$ position of each particle is determined through the iterative use of an accept-reject algorithm, accepted so long as the position of the particle maintains the desired surface density profile.The $z$ position of the particle is similarly determined, accepted so long as the position maintains hydrostatic equilibrium. The velocity of each particle is exactly Keplerian. This technique results in discs that are identical in their global properties, differing only through a small amount of noise at the particle-separation level.

All discs are evolved until it is no longer computationally feasible to continue, which in reality means the density of a fragment has become so high that timesteps cannot be computed without switching that mass to a sink particle. However, since here we wish to examine physical and orbital properties of the fragments which are influenced by their radial mass distribution, we do not do this. All of the simulations fragment to form at least two bound objects, and their ultimate configuration is shown in Figure 1, which shows 9 column density plots, in physical units, and illustrates a variety of fragmentation scenarios as the simulation’s final configuration. We discuss this in detail in section 3, but include images now to make the explanation of our methods in section 2.3 clearer.

2.3 Algorithms

We present here two methods of detecting fragments in SPH simulations, and one method of linking them together between timesteps. Once a fragment has been identified, we then refer to it as a "clump". The first method of detection is based on the clump finding approach of Smith et al. (2008) (which is in turn based on the publicly available CLUMPFIND algorithm developed by Williams et al. 1994). The basic idea is to perform an ordered search on SPH particles from high (physical) density to low density. The highest density particle $i$ forms the center of a clump, and if the next particle in the list is a neighbour (i.e. in close spatial proximity), it is also added to this clump. If it is not a neighbour, it forms the center of a new clump. This process is continued to the next most dense particle, until a minimum density threshold is reached.

The search in this manner, from least dense to most dense particle for our SPH simulations of protostellar discs was unsuccessful in identifying clumps in our simulation. We are faced with a different scenario to Smith et al. (2008), who used their algorithm in molecular cloud cores. Once our discs have evolved enough to fragment, the inner disc is so dense that many of the particles inside $\sim 20$ au fulfil the criteria to become the head of their own clumps. This results in neighbours, since this is a friends-of-friends algorithm, belonging to these even if they are in the outer disc, where fragmentation has actually occurred.

This problem was solved, to some extent, by using the gravitational potential of the particles, rather than the density, for the ordered search. We discuss this method in section 2.3.1. The inability of this method to identify low mass, fluffy fragments, or fragments that are so deep in the potential well of the central star that they are ultimately tidally destroyed, prompted the development of an approach that could correctly identify such fragments. The approach uses a gridded derivative search of the SPH interpolated density of the particles, and is discussed in section 2.3.2.

Finally, in both cases, the clumps are linked between timesteps using a merger tree algorithm, typical of dark matter halo tracers in cosmological simulations (see, for example, Srisawat et al. 2013). This process is detailed in section 2.3.3.

2.3.1 Gravitational potential search (CLUMPFIND)

Broadly speaking, clumps are created with a unique integer identifier (ID) at the local minima of the gravitational potential, so long as there are at least a minimum number of neighbour particles above some defined “noise” level. For our purposes, we define this critical number as $n_{\rm{critical}}=n_{\rm{mean}}-5n_{\sigma}$, where $n_{\rm{mean}}$ is the mean number of neighbours each particle has, and $n_{\rm{sigma}}$ is the standard deviation of the number of neighbour particles. We do this because the neighbour lists of particles at low density can be sparse due to the algorithms used to calculate the smoothing length, $h$.

We begin by creating a clump at the location of the central sink particle (star). All particles are then sorted by their gravitational potential energy, and we loop over the particles in order of most negative to least negative gravitational potential energy (i.e. most bound to least bound). We select the particle, $i$, with the most negative potential energy, and as long as it does not already belong to a clump (in which case, we exit and select the next particle), we iterate over particle $i$’s neighbours, $j$. If the majority of the neighbours $j$ (>50%) are in a clump $k$, the particle is also in clump $k$. We assign the particle ID of $i$ to ID${}_{i}=k$ and exit. If the majority of $i$’s neighbours are not in a clump, then since $i$ is the most bound particle it starts a new clump $l$, provided that $n_{\mathrm{neighbours}}>n_{\rm{critical}}$. We then loop over particle $i$’s neighbours $j$, assigning ID${}_{j}=l$ so long as $j$ is not already in a clump. We then proceed to next most bound particle, $i+1$, and repeat (Forgan et al., 2016).

2.3.2 Density Derivative Search

In this method, we compute a 2D grid in cylindrical $r$ and $\phi$ coordinates, and bin all particles into these grid cells. The maximum density of a particle in each cell is then taken to be the peak density in that cell, $\rho_{i}$, which gives us a 2D sheet describing density maxima. The number of bins required to reasonably identify all clumps varies due to the stochastic nature of the simulations, more fragments with low density at $\sim$ 100 au require a larger number of bins to properly resolve, with our resolution criterion being that the number of clumps ultimately detected by the search is equal to the number of clumps that are determined “by-eye”. If fewer clumps are detected by the search than “by-eye”, then the resolution is increased until we detect these low-mass, low-density clumps. A typical resolution is 10,000 radial bins and 7200 azimuthal bins.

We next take the derivative of the peak density in that cell with respect to $r$ and $\phi$. As this is noisy in density space, we smooth these derivative, equivalent to making a new signal, where the element is now the average of $n$ adjacent elements, such that:

[TABLE]

and

[TABLE]

Despite its simplicity, this approach is best at removing white noise while keeping the sharpest step response. The value of $n$ ($\sim$ 100 is typically sufficient), like the number of bins in which to bin the data, must be optimised by the user to get the best compromise between smooth data, which removes false peaks, and data which is sensitive enough to identify small clumps.

We then use this smoothed derivative to identify clumps, which will be a “real” peak in the density. A real peak is identified by a sustained zero-crossing of the derivative with a negative gradient. Peaks due to noise will also have zero-crossings, but they are sustained for fewer bins than real peaks. These false peaks can be eliminated by requiring that the zero-crossing is sustained for $m$ bins, with $m$ optimised by the user to remove most (if not all) false detections while still detecting less dense clumps.

The radial search is shown in Figure 2, which shows the radial density profile (blue solid line) and the derivative of the radial density profile (dashed red line) for the disc shown in simulation 6 of Figure 1. The zoomed section shows the peak of a clump, with the negative gradient zero-crossing of the derivative identifying the peak.

Once the particle marking the center of a clump ($i$) has been identified, we add all of $i$’s neighbour particles $j$ to that clump. We now loop over all particles which form that clump, adding their neighbours to this clump as well. We repeat this until we reach some density threshold. We found that adding neighbour particles until more than half the particles in the neighbour sphere are less dense than the inner 1 au of the disc produced good results. Once we have identified the bulk of the clump, we then proceed with a potential search described in section 2.3.1, which determines to which clumps the rest of the unidentified particles in the simulation belong.

2.3.3 Merger Tree

At this point, we have a set of clumps in each timestep, and we need to track them over the duration of the simulation. To do this, we use a standard algorithm from halo tracking in cosmological simulations (see, e.g., Springel et al. 2001). Each clump, in each timestep, is given an integer ID by our algorithm. So that we can trace the evolution of this clump throughout the simulation, these IDs must be linked. Since we are modelling fluid through the use of pseudo-particles, the particles, that make any given clump, change between dumps, sometimes substantially. To link clumps, the crucial factors are the most-bound particle MBP, and shared member fraction (SMF). In our density-derivative search, we actually trace the most-dense particle, rather than the most-bound particle, but we use the term interchangeably to avoid the introduction of unnecessary acronyms. To be identified as the same clump between timesteps, they must share the MBP and have an SMF of > 50%. In some particularly volatile simulations, when using the density derivative method, the MBP may change, and the SMF may be < 50%. In this case, some of the clumps need to be manually linked during post processing by tracing a group of particles in each clump in each timestep. The basic algorithm is as follows:

Loop over clumps $i$ in previous timestep lastdump. 2. 2.

Find clump $j$, in this timestep thisdump, which contains the MBP from clump $i$ in lastdump.

1. 2.i.

If MBPi does not belong to any clump in thisdump, clump $i$ is not present in thisdump. 2. 2.ii.

If MBPi belongs to clump $j$ in thisdump, and clump $i$ and $j$ share at least 50% of particles, then IDj = IDi, and the clumps are linked between the two timesteps. 3. 2.iii.

If MBPi belongs to clump $j$ in thisdump, but clumps $i$ and $j$ do not have SMF > 50% of particles, then clump $j$ keeps its present ID. 3. 3.

End loop over previous timestep. 4. 4.

Loop over clumps in thisdump, checking for clumps with no progenitors in lastdump. Increment the maximum number of distinct IDs by the number of new clumps, and assign each clump the correct ID.

3 Results

We ran a total of 9 SPH simulations with almost identical initial conditions, differing only in the random number seed used to initialise the disc. Column density plots of the 9 simulations are shown in Figure 1. Despite the almost identical initial conditions, there is a large variation in final configurations and number of clumps in the system.

Our results sections fall into three broad categories- first, the relative merits of the two methods and the difference in clump detection between them. We show that, interestingly, our density derivative search detects all clumps that are detectable by eye, whereas our ordered potential search does not. In fact, generally speaking, the potential search does not detect any clump that is eventually destroyed, giving a good indicator of the likelihood of a clump’s survival.

We next discuss the implications of our results for current population synthesis models, comparing our clump mass and semi-major axis functions to existing population synthesis models. We show the clump interaction needs to be included in GI population synthesis models as early as during the gas phase, as scattering plays an important role.

Finally, we discuss interesting events in the simulations themselves. We introduce a piece of nomenclature now, to avoid confusion, that S$a$C$b$ means simulation $a$, clump $b$. This abbreviation is given in the title of any plot of a specific clump. Note that our clump numbering begins at 2, since clump 1 is the star+disc system. We also state now, for clarity, that any mass stated for our clumps includes unbound material. This is deliberate, in order to track more of the mass of the clump. Furthermore, what is currently unbound material around the clump, at these very early times, may eventually lose spin angular momentum through interactions with material in the disc, ultimately becoming part of the clump. By including the unbound material, we trace more of this process from earlier times. For high mass clumps ($>20$ MJ), the amount of unbound material is small, typically around $10\%$ or so. This is larger for lower mass clumps, up to $\sim 25\%$ of material identified may be unbound, rising to 40$\%$ in clumps ocurring in particularly volatile simulations that have many clumps, as their formation is often disrupted by interactions with other.

We look at the orbital properties of the clumps, and discuss clump mergers and tidal destructions by the central star. We show that: (1) destruction and merging are fairly common, (2) interactions between clumps can result in a clump changing its direction of spin from prograde rotation to retrograde rotation, and (3) retro-rotating clumps typically have more dramatic changes in their radial temperature profiles than prograde clumps.

3.1 Relative merits of gravitational potential energy search and density derivative search

We identified clumps in our simulation using two different search methods, an ordered potential energy search (OPS) based on Smith et al. (2008), and a novel approach based on a 2D density derivative search (DDS). Our first conclusion is that searching using the ordered potential of the particles only detects clumps that survive the duration of the simulation. However, using the DDS method, destroyed clumps are also detected. This is shown in Figure 3, which shows the total change in semi major axis, from when the clump is first detected to when the clump is last detected (or destroyed), against the total time for which the clump exists, i.e. from initial identification to the end of the simulation, or the last timestep in which it is identified, if it is destroyed. Larger markers indicate more massive clumps. The left hand panel shows the DDS results, the right hand panel shows the OPS results. Circular markers indicate clumps that have survived until the end of the simulation, square markers indicate clumps that are destroyed, and triangle markers indicate a clump that merged into another clump. There are no identical markers in both plots because the different algorithms detect the clumps at different times, therefore they migrate different distances.

In addition to containing no destroyed clumps, the OPS sample also has a relative insensitivity to clumps which will have a final mass of less than $\sim$ 5 MJ, and detects most of the clumps later in the simulation. This is shown in Figure 4, which shows mass accretion histories for clumps in each of the 9 simulations, as is detected the OPS. By comparison, Figure 5 shows the same 9 simulations, but with the mass accretion histories of the clumps determined by using the DDS method. As can be seen, many low mass clumps evade detection entirely under the OPS method; for example, simulation 6 has an additional 3 clumps that are not detected by the OPS, and those that are detected are generally detected later, such as clump 2 in simulation 2, which is detected $\sim 400$ years later in the OPS than in the DDS.

This is due to the nature of the potential search. Figure 6 shows the radial gravitational potential energy profile of the disc in Figure 1, simulation 6. Since the OPS proceeds from the particle with the most negative gravitational energy to the most positive, the OPS detects the clump at $\sim$ 80 au in Figure 6 first, and then detects the clump at $\sim$ 375 au second. However, it fails to detect the clumps at $\sim$ 50 au, $\sim$ 125 au and $\sim$ 200 au. This is because particles at the potential energy of the main body of the disc are identified as belonging to the clumps with the most negative potential energy during the neighbour check described in section 2.3.1. Then, when the particles belonging to clumps at 50, 125 and 200 au are checked, they are found to already belong to either the main body of the disc (i.e. the central sink) or one of the clumps with the deepest potential well.

Although this could be fixed by adopting a gridded approach to the potential search (thereby eliminating the dominating effect from the clumps with the largest potential well), OPS has a desirable feature, namely demarcating clumps that are likely to survive the simulation, and those that are not.

The OPS method’s insensitivity to small clumps, and reliance on deep potential wells for identifying the body of the clump, mean that fewer clumps are detected by the OPS method, and those that are identified are often initially identified at artificially small masses ($\sim 10^{-3}$ MJ), as there is only a small amount of mass with a potential well deep enough to be identified. This can be seen in Figure 4, which shows mass accretion histories for all 9 simulations using the OPS method. In every simulation, what ultimately grows to be the largest clump is initially identified with a mass well below the Jeans mass. The DDS search method, however, does a better job of correctly identifying the mass associated with the young clumps, typically identifying between 5 MJ and 10 MJ of mass. This is shown in Figure 5, which shows the mass accretion histories for all 9 simulations, as identified using the DDS method.

In addition to the difference in measured initial clump mass, both methods differ in the final mass attributed to clumps. Typically, the growth is smoother for the OPS method, since once mass is deep in the potential well it is unlikely to change. However, if the gravitational potential energy profile of the disc changes then so too will the attributed mass of the clumps.

Another feature of the DDS method is its ability to identify low mass, low density clumps that do not have a strong signal in their potential. This can be seen in Figure 5, simulation 6, which identifies an additional three clumps compared to Figure 4, simulation 6. Comparing Figure 6 and Figure 2, we can see that these clumps have much stronger signals in their radial density profiles compared to their radial potential energy profiles. This is a useful predictive feature, since the left hand panel of Figure 3 and all of Figure 5 show that $\sim 20\%$ of the fragments in our simulations are tidally destroyed (we discuss the implications of this for population synthesis in section 3.2), none of which are detected in the OPS method.

Therefore, if a clump is detected in the DDS and not in the OPS, it is indicative that either the clump will stay relatively low mass and not accrete further, or that it will be tidally destroyed.

3.2 Comparison to gravitational instability population synthesis models

We ran the Forgan & Rice (2013) GI population synthesis (GIPS) models for 4000 years, which is comparable to the timescale for which our SPH simulations are able to run until no longer computationally feasible. We dictated that the opacity power law be $p_{k}=1$, and that the disc is not truncated after fragmentation. Figure 7 shows initial and final mass and semi-major axis distributions for the GI population synthesis model and for the clumps in our SPH simulation. The red, hatched histogram is the population synthesis, blue outline is the SPH clumps. The left hand panel of Figure 7 shows the initial and final semi-major axis distributions for both samples, and the right hand panel shows initial and final mass distribution for both samples. We cut off the tail of the initial and final masses beyond $>35$ MJ, since we cannot feasibly simulate masses above this without switching to sink particles.

Comparing the initial mass and semi-major axis distributions in the GIPS models to our SPH clumps, it would initially appear that clumps in SPH simulations form much further out, and at much lower masses, than in the population synthesis models. In reality, this is somewhat a limitation of our identification algorithms, as at very early times, the clumps can escape detection because of their low density/less negative gravitational potential energy, such that they have already undergone some radial migration before they are detected. If they have undergone sufficient radial migration, they may be far enough out in the disc to not accrete much material, hence remaining low mass.

Having established that the initial mass and semi-major axis distributions for our SPH clumps are subject to some limitations of our detection algorithm, we now compare the final mass and final semi-major axis distributions of our SPH clumps to those in the GI population synthesis model. First, the bottom left hand panel of Figure 7 shows a dearth of clumps at $R<25$ au, when compared to the GIPS model - this is simply due to the last measured value of $a$ before the fragment is destroyed. This figure also shows that the distribution of semi-major axes is different to what is expected, given the GIPS model data, and therefore the mechanism that allows these separations to exist at early times (i.e. clump-clump interactions) plays an important role in the ultimate orbital distribution function of the sample.

Second, the bottom right hand panel of Figure 7 shows that our final mass distribution is bimodal, with peaks at $\sim 5$ MJ and $\sim 30$ MJ. This is somewhat consistent with previous measurements of mass distributions of fragmenting discs by Vorobyov et al. (2013), who find that there are two maxima in their mass distribution, at $\sim 5$ MJ and $\sim$ 60 MJ. Unlike Vorobyov et al. (2013), we find a gap at $\sim$ 15 MJ, whereas they find a minima at $\sim$ 25 MJ, and our second maxima is at $\sim$ 30 MJ rather than $\sim 60$ MJ. Given our small $N$ statistics, we would probably expect our distribution to converge on a minima at $\sim$ 15 MJ, rather than the gap that is currently present. Additionally, our second peak is capped at $\sim$ 30 MJ in our simulations since this is typically when the density in a fragment becomes so high that it is computationally unfeasible to continue the simulation. With increased computation time, the mass of our largest clumps would probably increase.

Since our algorithm is quite robust at later times, detecting to within a factor of 2 the "by eye" clump mass in low-mass clumps (considered only in the bound region of the clump), we can see, in the final mass distribution shown in the bottom right panel of Figure 7, that the GI population synthesis model is significantly underestimating the fraction of planets at $<5$ MJ, even accounting for under-estimating low mass clumps by a factor of 2.

This can be explained by Figure 9, which shows the mass semi-major axis distribution for the final values of the SPH clumps and the population synthesis fragments. The SPH clumps are the dark circles, and the population synthesis fragments are the light squares. As can be seen, low mass clumps in our simulations are scattered out to large semi-major axis at these early times, and fragment-fragment interactions are likely to play an important role in the ultimate fate of a fragment. If it is scattered out to large $a$, it is much less likely to be tidally destroyed and far more likely to survive the duration of the simulation. This would suggest that GI population synthesis models need to include fragment-fragment interactions in this early gas phase, since current models suggest that $\sim 40\%$ of initial fragments are tidally destroyed. If a significant fraction of these are scattered out to large radii, their survival rate could potentially be much higher.

3.3 Orbital properties

We carry out analysis of the orbital properties of our clumps only using the sample as detected by the density derivative search, as this method is sensitive to most clump masses and semi-major axes. The total semi-major axis evolution of all clumps is shown in the left hand panel of Figure 3, which we have already discussed, and refer the reader back to. Circles mark surviving clumps (including clumps that subsume another clump), squares mark destroyed clumps, and triangles mark merged clumps. Larger markers correspond to more massive clumps. For destroyed clumps, we take the last measured mass. Roughly half of our most massive clumps migrate radially inwards, which is consistent with migration in locally isothermal discs, as objects exchange angular momentum with the surrounding gas and move inwards. However, about half of our most massive clumps migrate radially outwards. This is known to be possible in radiative discs (Kley & Nelson, 2012), but requires either large torques or steep surface density gradients (D’Angelo & Lubow, 2010). Large torques can have many sources, but in massive, self-gravitating discs they are likely to be in the form of global spiral arms. We carried out a Fourier analysis on the density structure of our discs, to determine the Fourier amplitude of each $m$-mode (where $m$ is the number of spiral arms). The amplitude, $A_{\mathrm{m}}$, of each mode, $m$, is calculated by

[TABLE]

where ${N_{\mathrm{region}}}$ is the number of particles in the region we are considering (for our case, $R=20$ to $R=100$ au), and $\phi_{i}$ is the azimuthal angle of the $i^{\mathrm{th}}$ particle. Some example amplitudes are shown in Figure 11, which shows the first 10 Fourier components of the density structure of 2 discs in their initial state (i.e. when they have just begun to fragment), marked in red, and the same 2 discs in their final state, marked in black. The discs are from simulation 1 and simulation 5, and their final state can be seen in their column density plots, shown in Figure 1. These discs were selected because they ran for the same length of time, and they have contrasting final $m$-modes states, simulation 1 ultimately peaks in the $m=2$ mode, and simulation 5 ultimately peaks in the $m=6$ mode.

In this fashion, we determine the dominant final $m$-mode in each disc, and plot the final semi-major axis of our clumps as a function of this $m$-mode, in Figure 9. We note that our decision to name a dominant mode, based on a relatively low amplitude difference, may be questioned. However, these discs are not in a quasi-steady state, having undergone fragmentation, so persistent spiral modes may be unlikely to form due to tidal disruption from these clumps. Despite the transient nature of the spiral modes, global, low $m$-mode spiral arms can exert considerable torque, and this is clearly important for the final orbital configuration of our clumps, as shown in Figure 9, which displays a rough empirical relationship between the maximum semi-major axis of a clump, $a_{\mathrm{max}}$, and $m$, such that

[TABLE]

Of course, this relationship is preliminary, since we only consider 9 discs, all of the same mass, and it has been shown that more massive discs are dominated by low $m$ spirals (Lodato & Rice, 2004, 2005). Indeed, since it has been shown that in discs without fragmentation we expect the number of spiral modes to be related to the disc-to-star mass ratio, $q$, such that $m\propto\nicefrac{{1}}{{q}}$ (Dong et al., 2015), to examine the full parameter space of spiral modes requires a range of $q$ values. Since we consider discs with identical $q$ values, why these discs have different dominant $m$-modes is a valid question. Again, this may be explained by these discs having fragmented into bound clumps. Bound objects in a gaseous disc produce stronger, more persistent torques than spiral density fluctuations alone. If one of our clumps is scattered out of the main body of the disc, say, by interaction with another clump, it may exert a tidal torque on the material in the disc (and by Newton’s third law, the material in the disc will also exert a force on the clump). Larger torques are associated with lower $m$-modes, and if these tidal responses from the discs to the clumps is responsible for the low $m$-mode domination in the disc, then we expect to see steeper surface density profiles in discs with low $m$-modes, since more mass will have been redistributed inwards as a result of this torque. This is difficult to unequivocally demonstrate in our set of simulations, since each simulation was run for a different length of time, and the amount of mass redistributed increases with time.

However, Figure 11 shows the mass enclosed as a function of radius for the two discs plotted in Figure 11, that ran for the same length of time. Although by no means conclusive, the slight increase in mass for a given radius between 40 au and 160 au for the $m=2$ disc in simulation 1 is consistent with tidal forces being responsible for the low $m$-mode becoming more dominant. To properly establish the nature of the preliminary relationship detailed in equation 5 therefore requires a range of disc masses and fragmentation scenarios, and we leave this to future work. Figure 12 shows the relationship between eccentricity, $e$, and semi-major axis for our SPH clumps. Larger markers indicate more massive clumps. For the most part, the more massive clumps are located on close in ($a\sim 50$ au), low eccentricity ($e\sim 0.1$) orbits, while lower mass clumps are at larger separations and with higher eccentricity. Since disc fragmentation forms objects on low eccentricity orbits ($e<0.1$), we can see the importance of clump-clump interactions in determining the final orbital properties of a clump. Very large eccentricities ($e\sim 0.7$) at large $a$ indicate that a clump is close to ejection, as excitations beyond unity ensure a clump is ejected from the disc. The top two panels of Figure 13 show the initial eccentricity distribution (left) and initial inclination distribution (right) of our SPH clumps. The bottom two panels of Figure 13 show the final eccentricity distribution (left) and final inclination distribution of the same. Inclination is calculated relative to the orbital plane of the central star, such that

[TABLE]

where $i$ is the orbital inclination of the clump, $\vec{L}$ is the orbital angular momentum vector of the clump (calculated relative to the centre of mass and centre of velocity of the central star), and $L_{\mathrm{z}}$ is the $z$ component of $\vec{L}$. Using least-squares regression, each plot has been fitted with a Gaussian of the form

[TABLE]

where $F_{\mathrm{p}}$ is the fraction of planets, $A$ is the amplitude of the curve (without offset), $x$ is either eccentricity or inclination, $\mu$ is the mean of the fitted distribution, $\sigma$ is the standard deviation of the fitted distribution and $k$ is the fitted offset constant. The fitted values are given in Table 1, and in the legend of each plot. The initial and final inclinations have been fitted with two distributions, the dotted blue line includes all points, and the solid red line does not include the most inclined point in each distribution, since there is a large gap between that point and the rest of the clumps, and a small sample size, it is unclear if this point is actually an outlier. We have provided these fits since current GI population synthesis models do not include orbital eccentricity or inclination, and, despite our small sample size, this information may be useful in the future development of these models. Aside from a small decrease in standard deviation, there is little change between our initial and final eccentricity distribution; both peak at $e\sim 0.1$ and share an offset constant of $k=0.006$. However, the orbital inclination of our clumps is reduced by a factor of $\sim 100$ between the initial and final states, showing that clump orbital inclination, in our SPH simulations, is rapidly reduced after formation. Considering that many of our clumps undergo dynamical interactions that cause scattering and eccentricity pumping on short timescales, this high degree of coplanarity may be surprising, especially when considering that most exoplanets have mutual inclinations of a few degrees (Figueira et al., 2012; Fang & Margot, 2013). However, it is consistent with our current understanding of highly inclined planet orbits relying on dynamical perturbations such as the Lidov-Kozai mechanism (Naoz et al., 2013). Our results may indicate that developing inclined orbits is difficult while a gas disc is present, even if substantial dynamical interactions between clumps take place in this time.

3.4 Spin properties

We analysed all of the fragments in our simulations, and found that several of them survive to the end of the simulation whilst undergoing retrograde rotation. This is shown in Figure 14, which shows the relative alignments between the orbital angular momentum vector and the rotational angular momentum vector of the clumps. Both the top and bottom panel is split into two parts, showing significant disalignment at the top, marked in red crosses, and good-alignment at the bottom. The top panel shows the alignment as a function of mass, and the bottom panel shows the alignment as a function of semi-major axis.

This prompts the question - how did they get to be retro-rotating? Did they form like this, or were they perturbed in some way? Having checked all of the clumps with significant disalignment, we can see that all of them were perturbed by a close encounter with another fragment, which typically flung them rapidly further out into the disc. We give the most extreme example in Figure 15, which shows the retrorotating clump 4 in simulation 4. In the leftmost panel, we see in the bottom left corner two clumps undergoing a close encounter. One of them then decreases its semi-major axis, whilst the other one is flung further out into the disc, to become retrorotating. Figure 17 shows the specific angular momentum profiles of a retro-rotating clump (top) and a prograde rotating clump (bottom). The blue line, at $T=2788$ years, is as soon is the clump is detected by our algorithm. We can see that the majority of the clump is in prograde rotation, with only the outer $\sim 1$ au in retrograde rotation. However, as time progresses and the clump continues to accrete material that is retro-rotating due to the encounter, this ultimately changes the rotation of the whole clump. For comparison, the bottom panel contains a prograde rotator of similar mass from simulation 5.

Figure 17 shows the rotation velocity profiles of two clumps. Top shows the ultimately retro-rotating clump 4 in simulation 4, and the bottom, for comparison, is a clump of comparable mass that is always undergoing prograde rotation. Negative velocity is determined by the orbital angular momentum vector of the whole clump and the rotational angular momentum vector of the material being anti-aligned. Both panels show the clumps at four different times, and the dashed lines indicate the breakup velocity profile of each clump. Both clumps are rotating under their breakup velocity for radii below 1 au, and their velocity profiles are consistent with solid body rotation (i.e. $v\propto R$) at these radii. Much further out, the clumps become nebulous, but we have included information out as far as possible to show the interesting angular momentum exchange between material at $T=2788$ years and $T=3377$ years for the retro-rotator. For the prograde rotator, the ultimate configuration is a good approximation to a solid body rotation curve out to extended radii ($\sim$ 3 au). Interestingly, for the retro-rotator, the velocity profile at the outer part of the fragment ($\sim 2.2$ au to $\sim 4$ au) is consistent with Keplerian rotation (i.e. $v\propto\nicefrac{{1}}{{\sqrt{R}}}$). This suggests the presence of a disc, or a disc-like structure, around the clump. Unfortunately, it is not (at the time of writing) currently possible to self-consistently re-resolve such regions in SPH simulations, so we are unable to investigate this further.

3.5 Density and temperature properties

Maximum density and temperature for all clumps in our simulation is shown in Figure 18, open circles show prograde rotating clumps and dark triangles show retrograde rotating clumps. The maximum temperature we identify in a clump is 2081 K, which implies that none of our fragments have started to dissociate molecular hydrogen and can therefore be considered as the first hydrostatic cores. Additionally, a further 9 clumps have internal maximum temperatures above $\sim 1000$ K, which means they would begin to evaporate dust. Both of these results are in good agreement with previous studies of fragments, such as Vorobyov et al. (2013), who found similar results. Our measurements of clump maximum temperature are limited by our resolution, since running simulations at higher resolution would allow higher densities to be reached before becoming computationally infeasible to continue.

As discussed in our introduction, the advantage of using a purely hydrodynamical simulation with no sink particles is that we can examine fragment internal structure during the simulation. With this in mind, we show the radial temperature and density profiles of 7 fragments, one of which ultimately becomes the hottest fragment in our simulations, shown in Figure 19, and the remaining 6 are three retro-rotators and three prograde rotators, of comparable mass and with similar simulation histories, shown in Figure 20. In both figures, each image is split into two panels, radial temperature profile on the top and radial density profile on the bottom. The initial clump profile is shown in solid red, and the final clump profile is shown in solid blue. For comparison, a polytrope is also plotted in each panel.

Figure 19 shows the radial density and temperature profile of the hottest clump identified out of our 9 simulations. We can see that although the radial temperature profile is well described by a polytrope of index $n=1.6$, for the initial state, and $n=1.8$ for the final state, in both cases the actual density deviates from the polytropic density significantly. This is also true for the rest of the clumps, shown in Figure 20, where the left-hand column shows prograde rotators, and the right-hand column shows retrograde rotators. In all cases in Figure 20, we have plotted a polytrope of index $n=1.5$, which is appropriate for a fully convective star such as a brown dwarf. Since in a polytrope, pressure $P$ and density $\rho$ are related by

[TABLE]

where $K$ is a constant and $n$ is the polytropic index, it implicitly assumes that pressure is a power law function of density which is constant throughout the star. For our clumps, it appears to be the case that a polytropic approximation may be too simplistic when estimating the internal structure of the clump. This may have some implications again for current GI population synthesis models, since orbital parameters are sensitive to the radial distribution of mass in the forming fragments. It is probably best to exercise caution when constructing clumps in GI population synthesis models, and it may be the case that a "follow the adiabats" approach, as used in Nayakshin & Fletcher (2015), is more appropriate.

When we compare the final states of the retro-rotating clumps with their prograde rotating counterparts, there is not much that would mark them as retro-rotating, final density and final temperature profiles for both directions of rotation are similar, and are consistent with other simulations of fragmenting protostellar discs (see, e.g. Vorobyov et al. 2013). There is one clump that is an exception when compared to the rest of the clumps, and is shown in the top right hand panel of Figure 20. This is S4C4. In its final state, the surface temperature of the clump is a factor of $\sim 3$ higher than we see in the rest of the clumps, and its central temperature is $\sim 425$ K, a factor of $\sim 2$ larger when compared to the rest of the clumps in Figure 20. This high temperature could possibly be explained by an encounter with another clump. Figure 15 shows this interaction, (clump 4 is marked by a green square). It is scattered by the more massive clump into the outer part of the disc, becoming shocked as it passes through a spiral arm. The encounter with the other clump causes the direction of rotation of clump 4 to change, but the large increase in temperature could be due to this motion through a region of increased density, entering perpendicular to the spiral arm where the density gradient is at its largest.

We suggest, therefore, that clump-clump interactions may provide a mechanism for dramatic increase in temperature of forming clumps, either directly through the interaction, or their subsequent scattering through dense regions in the disc. This mechanism could present a problem for terrestrial planet formation through tidal downsizing. For a terrestrial planet to form in the tidal downsizing scenario, a clump must not accrete too much mass, and then become tidally disrupted after migrating too close to the host star (Nayakshin, 2010; Boley et al., 2010). It then leaves behind a solid core, if dust grain sedimentation was sufficiently rapid to form a core. Since we have tentatively shown that clump-clump interactions are common (given our small sample size), then it may be possible that clump temperatures are frequently too high, at too young an age, for dust sedimentation to have taken place inside clumps. For a solid core to form, clumps need to exist at a temperature below the dust sublimation temperature ($\sim 1200$ K) for the duration of the sedimentation process.

3.6 Tidal disruption and mergers

Of the 41 clumps detected by the DDS method in our simulations, 7 were tidally destroyed by the central star, and 4 clumps underwent mergers. Despite our small sample size, these results suggest that both tidal destruction and mergers are common amongst protostellar disc fragments. An example merger is shown in Figure 21, which shows simulation 3, clumps 2 and 4, merging together. An example tidal disruption is shown in Figure 22, which shows simulation 7, clump 4 undergoing tidal disruption. Such tidal disruptions could potentially be an explanation for outburst type behaviour in young protostars, due to the rapid increase in accretion rate onto the central star (see, e.g.,Vorobyov & Basu 2005; Boley et al. 2010; Nayakshin & Lodato 2012). We leave observational signatures of our tidal disruptions to future work.

4 Discussion and conclusion

In this paper, we have presented one original method for identifying fragments in a simulation, the density derivative search, and one method adapted from the CLUMPFIND (Williams et al., 1994; Smith et al., 2008; Forgan et al., 2016) algorithm. We ran 9 SPH simulations of a $0.25$ M*⊙* disc around a 1 M*⊙* star, each with a surface density profile of $\Sigma\propto r^{-1}$, an inner radius of 10 au and an outer radius of 100 au. Each simulation was run for as long as computationally feasible without converting dense regions to sink particles, since we wished to calculate orbital properties for our fragments, which are sensitive to their radial mass distribution.

Each disc fragmented to form at least 2 bound objects, and we analysed the fragments (which we call clumps, once they have been detected) using the density derivative search and the adapted CLUMPFIND method. We have shown that these two methods are complementary, as the density derivative search is able to detect low mass clumps, and clumps that are ultimately tidally destroyed, while the search using the adapted CLUMPFIND method filters out clumps which are unlikely to survive the simulation, but also has a relative insensitivity to low mass clumps.

We compare our sample of fragments to the population synthesis model of Forgan & Rice (2013), and find that our algorithm has some limitations at early times (i.e., during the initial period of fragment formation), which means that some radial migration has already happened before the algorithm detects the fragment. Despite this, it is fairly robust at late times (i.e., a few orbital periods after formation), and so our final mass and final semi-major axis functions, and the mass semi-major axis relationship, are representative of the final configurations (i.e., at $T\sim$ 4000 years) of our systems. Of course, the final mass and final semi-major axis distributions that we present here will differ from the ultimate distributions of the systems. These will only be determined some $\sim 10^{6}$ years after formation, long after the disc has dispersed. What happens after disc dispersal is not considered in this work, but we refer the interested reader to Forgan et al. (2015) and Li et al. (2016).

We examine the internal temperature and density structure of our fragments, and compare them with appropriately indexed polytropes. We find that the central density of our fragments are typically an order of magnitude denser than their polytropic equivalent, and since orbital parameters of a body are sensitive to its internal mass distribution, we recommend caution when using polytropes in GI population synthesis models to calculate orbital parameters.

Furthermore, that the interiors of the clumps may not be well described by polytropes may raise concerns about the validity of the Forgan et al. (2009) hybrid radiative transfer method, which uses the polytropic cooling formalism of Stamatellos et al. (2007) to approximate radiative cooling. Essentially, the assumption is that each SPH particle is embedded in its own polytropic pseudo-cloud, of polytropic index $n$. Therefore, the cooling of any morphology that is not a polytrope of index $n$ may not be exactly described by the polytropic cooling formalism. However, so long as the geometry is approximately spherical, which gravitationally bound clumps tend to be, then the polytropic cooling approximation is valid to within an order of unity for a variety of opacity laws (Wilkins & Clarke, 2012).

More generally, the polytropic cooling method of Stamatellos et al. (2007) has been criticised when employed in planar geometry, such as in self-gravitating protostellar discs. The most alarming of these is the systematic under-estimation of cooling rates by a factor of $\sim$ 100 (Wilkins & Clarke, 2012), due to the use of local variables only in optical depth determination. This ignores the decreasing column density normal to the midplane of the disc, which may offer an easier escape path for photons.

However, this criticism is specific to the polytropic cooling method of Stamatellos et al. (2007), rather than the hybrid method of Forgan et al. (2009). By including flux limited diffusion, the problems outlined in Wilkins & Clarke (2012) are addressed.

More recently, the polytropic cooling approximation of Stamatellos et al. (2007) has been updated to include radiative feedback from sink particles (Stamatellos, 2015), which is not captured by the Forgan et al. (2009) hybrid scheme. However, we do not use sink particles in this work, since we aim to characterise the radial properties of our fragments. By using only SPH particles, the flux limited diffuser in the Forgan et al. (2009) algorithm naturally provides radiative feedback from the fragments, as heat is diffused from hotter to cooler regions.

Ultimately, the hybrid method of Forgan et al. (2009) is an approximation to a computationally much more expensive task, a full radiation hydrodynamical treatment of disc evolution. Therefore, results obtained using this approximation should be adopted with knowledge of its limitations in mind.

The radiative cooling of real protoplanetary discs depends strongly on opacity, which is not examined in this work. Dust opacities have a large degree of uncertainty; for temperatures below $\sim 1000$ K, dust grains dominate absorption, so this uncertainty is likely to matter in the region of parameter space considered in this work. For example, Semenov et al. (2003) have shown that for temperatures $100\,{\mathrm{K}}<T<1000{\,\,\mathrm{K}}$, opacity may differ by a factor of $\sim$ 10. Since it is known that initial fragment mass depends on opacity (Masunaga & Inutsuka, 1999), our initial fragment mass distribution may differ from what we state here.

However, it is widely accepted that dust opacity is proportional to metallicity, and it has been shown that while low-metallicity fragments may be a factor of $\sim$ 3 more massive than solar metallicity fragments, the initial mass distribution of the fragments is, apart from the shift at low mass, very similar between the two metallicities (Bate, 2005). We conclude, in light of this, that our fragment mass distributions are probably reasonable, but caution that the lower mass fragments may be larger if opacity is decreased. It is possible to estimate the mass dependence on opacity by appealing to the opacity limit for fragmentation, as follows.

Usual derivations of the opacity limit for fragmentation include an efficicency factor, $e$, which subsumes, amongst other unknowns, the opacity of the gas. Adjusting this efficiency factor allows us to examine a range of potential fragment masses. We begin with an expression for the power of a ball of gas collapsing in freefall:

[TABLE]

where $B$ is the gravitational binding energy of the gas, given by

[TABLE]

and $t_{\mathrm{ff}}$ is the free-fall time of the gas, given by

[TABLE]

For the collapse to continue isothermally, then

[TABLE]

where $L$ is the luminosity describing the radiation of the ball of gas, as a blackbody, given by

[TABLE]

where $T$ is temperature, $\sigma$ is the Stefan-Boltzmann constant, and $e$ is the efficiency factor of the radiation. Substituting our expression for power, $|B|/t_{\mathrm{ff}}$, and luminosity, $L$ into equation 12, we can rearrange equation 12 to arrive at an expression for the critical mass:

[TABLE]

above which the collapse proceeds adiabatically. However, this mass must be above the local Jeans mass, $M_{\mathrm{J}}$, for the collapse to continue, so we now have a condition

[TABLE]

which must be satisfied. If $M_{\mathrm{J}}$ and $M_{\mathrm{crit}}$ are equal, we then arrive at the opacity limit for fragmentation, which describes the minimum mass of a fragment that may form. Inserting constants, we now have an expression for the minimum fragment mass,

[TABLE]

From this, we can see that for a fixed temperature, minimum fragment mass increases with a decrease in efficiency. If the efficiency of the radiation is 100%, i.e., the gas has a sufficiently low opacity such that it is effectively completely optically thin to escaping radiation, then for a temperature of $T\sim 10$ K a fragment mass of $\sim 10$ MJ is expected. For an efficiency of 20%, then this rises to $\sim 23$ MJ. Clearly, how strongly this efficiency depends on opacity will determine the relationship of proportionality, however, it is clear that opacity plays an important, if not dominant, role in determining initial fragment mass. In addition to determining initial fragment mass, the opacity will also play an important role in determining the final fragment mass through accretion rate onto the fragments, since the thermodynamics of the gas inside the Hill sphere of the fragment depends on this opacity (Stamatellos, 2015).

We find that fragment-fragment interactions play a substantial role in the ultimate fate of our systems. Low-mass fragments can be scattered out to large radii, (and therefore remain low mass, since there is less material to accrete), and are therefore unlikely to be tidally destroyed by the central star. Since current GI population synthesis models suggest that $\sim 40\%$ of initial fragments are ultimately tidally destroyed (Forgan & Rice, 2013), if a significant number of these fragments are actually scattered out to large separations by interactions with other fragments, this figure could potentially be much lower. We therefore recommend that fragment-fragment interactions in the gas phase of the disc be included in any new GI population synthesis models.

During their lifetime, we find that our fragments may be shocked as they pass through spiral arms, rapidly increasing the internal energy of the fragment. This could have implications for terrestrial planet formation through the tidal downsizing hypothesis, since solid core formation requires rapid dust sedimentation. If the interior of these fragments are hot enough to sublimate dust at very early stages in their lifetime, it may not be possible for them to form solid cores, a key assumption in the core-assisted gas-capture hypothesis of Nayakshin et al. (2014).

In addition to encounters scattering fragments out to large radii, we also find a tentative relationship between the dominant azimuthal wavenumber in the disc, and the maximum semi-major axis of a clump in that disc, such that $a_{\mathrm{max}}\propto\nicefrac{{1}}{{m}}$. This seems to suggest that the spiral arms of protoplanetary discs play as large a role in the dynamical fate of the clumps as do the clump-clump scatterings. Although this relationship is preliminary, and requires more simulations in a range of disc-to-star mass ratios to confirm it (since $m\sim\nicefrac{{1}}{{q}}$ (Dong et al., 2015), where $q$ is disc-to-star mass ratio), it is unsurprising that such a relationship may exist, and we suspect that the relationship is one of inverse proportionality, rather than one of negative proportionality, due to the relationship between gravitational torque $\Gamma_{\mathrm{G}}$ and azimuthal wavenumber as follows (for a full derivation and comprehensive explanation, we refer the reader to Binney & Tremaine 2008). The gravitational torque $\Gamma_{\mathrm{G}}$ exerted on material outside a radius $R_{0}$ in a disk is given by

[TABLE]

where $\Sigma_{1}$ is a gentle function of radius, $k$ is the radial wavenumber defined as

[TABLE]

$mf(R)$ is the radial shape function of the spiral (for a simple example, see Hall et al. 2016), and $\mathrm{sgn}(k)=+1$ for trailing spirals (i.e, positive torque exerted outwards). At a given value of $R_{0}$, and the same function for $\Sigma_{1}$, we can see that

[TABLE]

To establish the relationship $a_{\mathrm{max}}\propto\nicefrac{{1}}{{m}}$, we have assumed that the amount of torque is directly proportional to the change in radial distance. Whilst this is almost certainly an over-simplification of matters, we can see from

[TABLE]

that if we assume mass, velocity, and the angle between $\vec{v}$ and $\vec{r}$ stay fairly constant, then $\Gamma_{\mathrm{G}}\propto\nicefrac{{dr}}{{dt}}$, i.e. the distance we wish to move our fragment. Assuming that all fragments form at roughly the same $r$ in a given disc, (i.e. where that particular disc becomes susceptible to fragmentation), then we recover $a_{\mathrm{max}}\propto\nicefrac{{1}}{{m}}$.

Varying the surface density profile in a protostellar disc will alter disc torques and migration rates of planets (see, for example,Baruteau et al. 2011). Since we simulated 9 protostellar discs that had initially identical surface density profiles, this is not something we have investigated, but it could potentially alter our results. For example, steeper surface density profiles correspond to more rapid migration rates. It is therefore feasible that a sharp cut-off at the disc outer edge, as we use in our initial conditions, may exacerbate migration torques in the disc. However, since the majority of the disc mass is contained well within the disc’s numerical outer edge, it is probably reasonable to assume it is not the dominant effect when considering the radial migration of fragments.

Despite the relatively short timescales of our simulations (since we did not make use of sink particles), we can see that orbital properties of fragments are drastically altered by interactions with each other. Since disc fragmentation forms objects with initially low eccentricities ($e<0.1$), it is generally accepted that measurements of eccentricity as a function of orbital distance will constrain the formation mechanism of giant planets and brown dwarfs (Vorobyov, 2013), with high eccentricity being caused by dynamical scattering.

However, our results in Figure 12 possibly suggest that these high eccentricity orbits could be formed at very early times, during the gas phase of the disc, and as such eccentricity measurements of brown dwarfs and giant planets may not, necessarily, constrain their formation mechanism. We have shown that the intial orbital inclination of our fragments is reduced by a factor of $\sim$ 100 over the duration of the simulation (Figure 13), despite the significant dynamical interactions many of the fragments experience. This suggests that although dynamical interactions certainly can create highly inclined orbits, doing so while the gas disc is present may be much more difficult.

On the other hand, our simulations are of discs in isolation. Inclination and eccentricity may be excited by environment, such as a stellar companion, or location within a cluster environment (Forgan et al., 2015). Since current GI population synthesis models do not include eccentricity or inclination information, we have provided several Gaussian fits in Figure 13 from our SPH simulation data. Despite our small sample size, we hope these plots will be useful in further development of GI population synthesis models.

In this work, we do not consider how solid particles will effect the formation, evolution, and survival of the clumps, nor their effects on the behaviour of the gas. This is, quite possibly, the most important limitation of our work. As we have already discussed, opacity could heavily influence the initial mass distribution of our clumps. It is also feasible that altering the opacity in the clumps, due to movement of solids, could result in clumps cooling more (or less) rapidly, resulting in more (or fewer) clumps surviving (Nayakshin, 2010). More clumps surviving because of local opacity changes could present a partial solution to the rapid inward migration and subsequent disruption of GI fragments found in so many other works (see, for example, Vorobyov & Basu 2005; Baruteau et al. 2011; Zhu et al. 2012).

To date, no global 3D SPH simulations of dusty, self-gravitating, fragmenting protostellar discs have been conducted. This is a computationally expensive task, with complicated physics. Almost certainly, a full treatment of radiation hydrodynamics would be necessary to correctly capture the cooling of the system, and the spatial evolution of the dust within the gas. High resolution is needed in the core of these fragments to trace their compositional change. However, if we are considering the global evolution of the system, we must be able to capture the full dynamical range. This is probably not possible with current computational abilities, and is one of many multiphysics problems present in protostellar disc modelling (Haworth et al., 2016).

How fragmentation proceeds when solid grains are included in the simulation is clearly an interesting question. Although not considered in 3D hydrodynamics simulations, this idea has recently been explored by Nayakshin et al. (2014), in their so-called "core-assisted gas capture" (CAGC) paradigm. In this scenario, grain sedimentation inside a fragment forms a core of heavy elements. Upon reaching a critical core mass, the surrounding gaseous envelope collapses onto the core, analogous to the core accretion paradigm. The assumption in the work of CAGC is, of course, that these cores do form. Gravitational instability tends to occur beyond $\sim 50$ au in protostellar discs (Rafikov, 2005). For a significant core to form in GI clumps due to grain sedimentation, there must be a substantial local dust-to-gas ratio at the site of clump formation. Since dust grains tend to migrate rapidly inwards (Weidenschilling, 1977), it may well be difficult to maintain a significant dust-to-gas ratio at the site of fragmentation, making core formation from grain sedimentation inside a fragment difficult. However, GI is also a rapid process, so as long as the fragmentation takes place rapidly, then the local change in pressure gradient may be sufficient to prevent inward migration of these grains.

If massive cores in GI clumps do form, as assumed by CAGC, then we should probably expect differences in the final clump distributions to what we have presented here. An interesting feature of the (CAGC) paradigm is that if planets are formed by this method, then we should expect a positive metallicity correlation.

Assuming, for now, that the collapse proceeds with no core mediation, then it is probably reasonable to expect no metallicity corellation; indeed, there is at least some numerical evidence that we may expect a negative metallicity corellation (Cai et al., 2005), since low metallicity corresponds to faster cooling (and therefore stronger spiral amplitudes, increasing the effective gravitational stress) which may result in fragmentation. This is because higher metallicity results in better cooling only in the optically thin regime, for example, at the tenuous surface of the disc.

Conversely, higher metallicity results in slower cooling in the optically thick regime, i.e. at the disc midplane, which is the location of fragmentation. Therefore, gas-giant planets formed through GI should be preferentially found around low metallicity stars. It seems to now be fairly clear that gravitational instability rarely forms planetary mass objects (Rice et al., 2015), but if it does, these objects will be much larger than Jupiter. We should, then, find that planets more massive than Jupiter are more frequent around low-metal stars. However, it has been suggested that this is not the case, since planets more massive than $\sim 3$ MJ seem to be found less often around low-metallicity stars (Thorngren et al., 2016).

Overall, the picture is unclear. While CAGC may result in a positive metallicity correlation, there seems to be enough evidence to suggest that fragmentation should preferentially occur in discs around low-metallicity stars, which, in turn, suggests that GI planets would be preferentially found around low-metallicity stars. In the CAGC paradigm, fragmentation must first take place, before the core can form through subsequent grain sedimentation. If this fragmentation happens less often around high-metallicity stars, this could still result in an over negative metallicity corellation. Ultimately, the ambiguity of what happens at low-metallicity does not remove the strong high-metallicity correlation with planet frequency. This is mostly suspected to be the work of core accretion, and is difficult to explain with any GI theory alone.

The work we have presented here outlines the first direct comparison between a GI population synthesis model and a suite of 3D, global hydrodynamics simulations of fragmenting protostellar discs. While the work of Forgan & Rice (2013) presents the first attempt at a GI population synthesis model, it is not the only one (see, for example, Galvagni & Mayer 2014; Nayakshin & Fletcher 2015, and a recent extension to the Forgan & Rice 2013 model in Forgan et al. 2015). As mentioned in the introduction to this work, our simulations are run for as long as is computationally feasible without the use of sink particles, since part of our aim was to characterise the internal density and temperature profile of these fragments. However, doing so limits us to very early times in the disc, typically around $\sim 4000$ years or so after the disc has initially formed. For this reason, we did not compare our results to the models of Forgan et al. (2015), which consider the ultimate dynamical fate of the fragments after disc dispersal.

Similarly, the work of Forgan & Rice (2013) considers the fragmentation phase of a disc, unlike Galvagni & Mayer (2014). Since our hydrodynamical simulations are analysed around the fragmentation phase, we wished to use a model that did not assume already that clumps exist in the disc, ruling out the Galvagni & Mayer (2014) models. Finally, we considered the Nayakshin & Fletcher (2015) models unsuitable for direct comparison to our hydrodynamics simulations because only one fragment per disc is simulated, unlike the population synthesis models of Forgan & Rice (2013), which places multiple fragments in a disc susceptible to fragmentation at separations of a few Hill radii. We stress here that we are not suggesting the superiority of the Forgan & Rice (2013) models, simply that those particular models were best suited for direct comparison to our hydrodynamics simulations.

Of the 41 clumps that are detected in these simulations, 7 were tidally destroyed ($\sim 20\%$), and 2 have orbits with eccentricity approaching unity ($e\sim 0.75$), which suggests that they are on their way to being ejected ($\sim 5\%$). If these clumps are ultimately ejected, then gravitational instability could, perhaps, also contribute to the population of free-floating planets (Rice et al., 2003b; Forgan et al., 2015). We have demonstrated that the orbital and structural evolution of neighbouring fragments are linked; we recommend, therefore, that any future population synthesis models are able to account for this.

5 Acknowledgements

We thank the anonymous referee for their constructive comments, which greatly improved the discussion in this paper. We would like to thank Daniel Price for his publicly available SPH plotting code SPLASH (Price, 2007), which we used to produce Figures 1, 15, 21 and 22. We also thank Fabian Fischer for his online LaTex Overlay Generator (available here: https://ff.cx/latex-overlay/), which we used to annotate Figures 15, 21 and 22. KR gratefully acknowledges support from STFC grant ST/M001229/1. DF gratefully acknowledges support from the ECOGAL project, grant agreement 291227, funded by the European Research Council under ERC-2011-ADG. The research leading to these results also received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement number 313014 (ETAEARTH). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 681601).