Formation of satellites in circumplanetary discs generated by disc instability

We investigated the formation and evolution of satellite systems in a cold, extended circumplanetary disc around a 10 $M_{\rm{Jupiter}}$ gas giant which was generated by gravitational instability. We use a population synthesis approach, in which we seed the disc with satellite embryos and let them migrate, accrete mass, collide, with subsequent generations being created until the disc is dissipated. In each run we choose randomly the dust-to-gas ratio, dispersion- and refilling time-scales, the number of embryos and their starting locations in a plausible range. The disc structure is the result of a 3D global disc SPH simulation \citep{Szulagyi16b}. We also investigate the effect of the planet's semi-major axis on the resulting satellite systems, taking 50 AU as the nominal case. In the nominal case we find that most satellites are close in mass to the Galilean ones, with a maximum at around 3 $M_{\rm{Earth}}$ and form typically on time-scales comparable to the dispersion time-scale. We also find that about 10 $M_{\rm{Earth}}$ worth of satellites migrate into the planet, polluting it with metals. The influence of a different planet semi-major axis is mainly felt through the change in disc size. We find that for the discs closer to the star, the satellites are lighter, formation time-scales are longer, and more moons are lost into the planet, as a smaller disc means easier migration for the massive satellites. The probability of detecting satellites like the ones in this work is very low ($\leq$ 3%) even with upcoming powerful telescopes like E-ELT.


INTRODUCTION
There are two major scenarios of gas giant planet formation, core-accretion (CA; Pollack et al. 1996) and gravitational instability (GI; Boss 1997,Mayer et al. 2002. In both formation mechanisms a circumplanetary disc (CPD; e.g. Lubow et al. 1999, Kley 1999, Shabram & Boley 2013) is created at the last stage of the formation process. Furthermore, in a similar way as planets form within the circumstellar disc, satellites can form in these CPDs (Lunine & Stevenson 1982). So far, a comparative study of the properties of the satellite systems formed in CPDs arising in CA or GI has never been attempted.
In our Solar System, every gas-and ice giant has a satellite system. However, we only have one possible exomoon candidate around an exoplanet so far (Teachey & Kipping 2018). Of course, detecting exomoons is particularly challenging, not only because of their tiny mass and size. It is E-mail: cassandra.inderbitzi@uzh.ch also because a priori, CPD forming massive planets are created beyond the snowline and these are the ones which can have satellite systems. Most exoplanet/exomoon detection methods are sensitive to planets much closer in. While dynamical interactions could bring giant planets closer to the star, during these some/all moons could detach from the planet or be torn apart in tidal interactions. Therefore, detecting exomoons is particularly challenging as the bulk of the population is orbiting in the outer planetary systems.
To understand the possible differences between satellite systems that form in GI or CA CPDs, the characteristics of these discs have to be taken into account. In the GI case, previous works have found massive CPDs: Galvagni et al. (2012) with disc masses of 0.5 M planet , while Shabram & Boley (2013) reported 0.25 M planet . However, given that the CPD is fed by the circumstellar disc through the meridional circulation (Szulágyi et al. 2014;Fung & Chiang 2016), it is natural that the CPD mass will scale with the CSD mass. On top of that, for GI to operate one needs a very massive CSD, which further contributes to the massive CPDs.
It was shown by Szulágyi (2017) that the CPD mass indeed linearly scales with the CSD mass, and it also scales with the planet mass. Regarding the temperatures of CPDs, the GI scenario predicts very low values: ∼100 K or below (Galvagni et al. 2012;Shabram & Boley 2013). CPDs that formed during the CA scenario were studied more frequently. The subdisc masses were spanning a large range from 10 −4 to 10 −2 planet masses (D'Angelo et al. 2003;Gressel et al. 2013;Szulágyi et al. 2014Szulágyi et al. , 2016bSzulágyi 2017), since again, the CPD mass will depend on the planet mass and the CSD mass. The temperatures, in works employing radiative simulations, were significantly higher than those of CPDs in GI simulations. Furthermore, within the CPD there was a large temperature gradient, ranging from several thousand Kelvin near the planet, to a few hundred Kelvin at the edge of the CPD (Ayliffe & Bate 2009;Szulágyi et al. 2016a;Szulágyi 2017;Szulágyi & Mordasini 2017). To understand the real differences between GI and CA CPDs, Szulágyi et al. (2016b) compared the masses and temperatures within two simulations for the two cases. The initial parameters of these simulations were the same, the planets were of 10 MJupiter at 50 AU from their stars. While the CA CPD was only 8 times smaller in mass than the GI CPD, their temperatures were an order of magnitude different, with the CA CPD being significantly hotter. This naturally rises the question of how different would be the satellite system formed from such markedly different conditions. So far, satellite formation was mainly studied by population synthesis. Sasaki et al. (2010) has investigated the outcome of Saturnian and Jovian moon systems based on the differences in disc mass, cavity between the disc and planet, as well as the termination time of gas infall. Miguel & Ida (2016) used a minimum-mass sub-nebula model for the CPD to reproduce the Galilean satellites. Their results have showed that the presence and location of the snowline could explain the different ice-rock compositions of the Jovian moons. Fujii et al. (2017) population synthesis very successfully reproduced the Galilean moons' 1:2:4 mean resonance. Instead of an assumed disc profile, Cilibrasi et al. (2018) used a hydrodynamical simulation on Jupiter's formation to serve as an input for the CPD's density-and temperature profiles. In this case, the authors also included mass influx from the circumstellar disc (Szulágyi et al. 2014) and a realistic dust coagulation model was used to determine the dust density profile (Dr ażkowska , where the satellite seeds were forming via streaming instability. While Shibaike et al. (2017) pointed out that the dust drift is an issue to keep solids and create satellitesimals, Dr ażkowska &  showed that dust traps and the outflowing gas stream can mitigate this effect and therefore satellitesimal formation is not only possible, but it happens very efficiently due to the short orbital time-scales around the planet. Szulágyi et al. (2018) used hydrodynamic simulations of CPDs of Uranus and Neptune, and population synthesis for their satellite formation to understand how those satellite systems could have formed. They found that the CPD that naturally forms during the formation of these ice giants can form Uranian-like moons systems. In the case of Neptune, such systems might have existed before Triton capture happened and wiped away the original satellite system. Recently, Shibaike et al. (2019) suggested that pebble accretion can reproduce many characteristics of the Galilean satellite system. Ronnet & Johansen (2019) successfully explained the formation of satellite systems around Jupiter and Saturn via pebble accretion as well.
Given the prospects of upcoming exomoon discoveries, here we study satellite formation with a population synthesis approach in GI generated CPDs from Szulágyi et al. (2016b). Given that the CPD characteristics, such as the disc radius as well as the temperature and mass influx from the circumstellar disc changes with semi-major axis of the planet, we rescaled our nominal CPD simulation at various orbital separations from the star. We compare the satellite systems with CA models as well as giving detectability predictions of the formed exomoons.

METHODS
In this work we perform population synthesis (Section 2.2), on a base of hydrodynamical simulation (Section 2.1) obtained disc profiles.

Hydrodynamical simulations
For the circumplanetary disc density-and temperature profiles, we use a hydrodynamical simulation output from Szulágyi et al. (2016b). This is a 3D global disc Smoothed Particle Hydrodynamic (SPH) simulation to study planet formation via disc instability. It is designed to mimic the formation of the HR 8799 system, which has several gas giants on wide (R > 30 AU) orbits, that might have formed via gravitational instability.
The circumstellar disc radius initially is set to 5-200 AU and 0.6 M Solar mass. The initial surface density profile is set up in such a way that it follows a power law, where the exponent is close to −1 in the region of about 30 -100 AU. For the temperature profile, an iterative procedure is used, which includes full force balance and stellar irradiation at t = 0 as well as disc self-gravity Rogers & Wadsley 2011). The stellar mass is 1.35 M Solar . There are several clumps that form within the circumstellar disc. We extract one of the smallest clumps, where the planet mass is ∼ 10 MJupiter with an equal mass CPD.
The simulation is run with the ChaNGa Tree+SPH code and contains 42 million particles. The finest resolution is 0.01 AU. Monaghan viscosity is used with α = 1 and β = 2 . The code uses radiative cooling, which is dependent on local gas properties, and gives an energy loss per time per volume like: where σ is the Stefan-Boltzmann constant, s = (m/ρ) 1/3 , Tmin is the minimum gas temperature background (around 10 K in this case) and τ the optical depth. This way, cooling is most efficient at around τ = 1, while it also recovers the dependence on τ of the cooling rate in the asymptotic limits of optically thin and thick discs. It does not solve the radiation hydrodynamic equation and as such does not take into account the accretional luminosity of contracting clumps, but it includes compressional heating (which is generated by PdV work) and shock heating.
Optical depths are computed using the tabulated values from D' Alessio et al. (1997) for the Rosseland mean opacity and from D' Alessio et al. (2001) for the Planck opacities, assuming solar metallicity and dust-to-gas ratio 0.01. To take into account the ortho/para ratio of molecular hydrogen, which varies as a function on temperature and is important to recover the thermodynamics across spiral shocks that occur in self-gravitating unstable discs (Podolak et al. 2011), a variable adiabatic index is used.

Population Synthesis
In this 1D semi-analytical model, the initial disc profile is taken from the hydrodynamic simulation of Szulágyi et al. (2016b).
Within the disc, satellite seeds are created, which can migrate, accrete, get captured into resonances, collide and get lost into the planet. At the same time, the disc is evolved too, both in density and temperature. A fixed grid is used within the CPD, with a cell size of a Jupiter radius and a fixed timestep dt = t disp /1000, where the dispersion timescale t disp is a typical time-scale on which the disc disperses. This parameter is set at the start of the simulation (see Section 2.2.3). The simulation ends once the system has reached a steady state, which is always before 14 t disp .

Disc structure and evolution
Because the population synthesis is performed in 1D, we azimuthally average and vertically integrate the CPD from the 3D hydro simulation. The initial density profile is shown on Fig. 2, while the initial temperature profile is on Fig.  3. Fig. 1 shows a density map of the CPD midplane and a vertical slice.
The CPD ranges from 3 RJupiter (as the planet is slightly bigger than 2 RJupiter and the grid size is taken to be 1 RJupiter) to 12555 RJupiter, i.e about 60 % of the planet's Hill-radius or 6 AU (according to the hydro simulation, see Section 2.1), and has a viscosity α = 0.004, equal to the value in the hydrodynamic simulation. The code units are RJupiter for the radius, MJupiter for the mass, Kelvin for the temperature, and years for the time. The initial profiles for gas density and temperature are fits to the simulation data and take the following form Σgas,0(r) = 4.99 · 10 −7 · e −2.32· r 1 AU The dust density is assumed to have the same shape as the gas density, but multiplied by a constant dust-to-gas ratio, which we change in the different 1D runs.
The temperature, gas and dust densities also evolve with time. For both, an exponential decay (Ida & Lin 2008) is assumed on a time-scale t disp : Where Tmin is the temperature of the gas surrounding the CPD, with a value of 11 K.
While this CPD is relatively massive, it is not selfgravitating, because the Toomre Q parameter (Toomre 1964) is always above unity. This is consistent with the absence of spiral structure and the large thickness of the CPD ( Figure 1).

Satellite formation and evolution
2.2.2.1 Satellite formation In this model a similar approach is used as in Miguel & Ida (2016), which means that seeds are created in the disc with a mass of 10 −7 MJupiter. The first generation is created between [0, t disp /2], and on positions between 1% − −80% of disc radius. They are uniformly distributed between [log10(0.01 · R disc ), log10(0.8 · R disc )]. A log-uniform distribution is chosen because the disc radius spans multiple orders of magnitude.
The model also allows for satellites to be created at later times, in a sequential manner, when one satellite is lost into the planet (Canup & Ward 2002). Using the results from Shibaike et al. (2017) to estimate the time-scale on which these later generation seeds form, a growth-time of the satellitesimal is assigned randomly between [10 4 , 10 5 ]. This satellitesimal growth is also exponentially suppressed on the dispersion time-scale, as on such a time-scale the dust reservoir from which satellites accrete, being proportional to the disc mass, will also undergo an exponential reduction.

Migration
Both Type I and Type II migration are considered and the two the regimes are separated by gap opening. At first, when the satellites are not too massive, Type I migration will occur. As we assume circular and coplanar orbits, migration manifests as a change in the orbital radius a and the migration prescription given by Tanaka et. al. (2002) is used: where h = cs Ω k is the scale-height of the disc, ΩK is the Keplerian velocity at radius a, cs is the local speed of sound, Msat is the satellite mass and M planet is the planet mass. The factor bI is a correction factor that depends on the disc for which the prescription by Paardekooper et al. (2010) is used in this work. The influence of the satellite on the gas density around it is also taken into consideration with a partial gap model. This is done by multiplying the velocity with a factor between 0 and 1, representing this depth, based on an analytical model by Crida & Morbidelli (2007).
As satellites increase in mass, they might entirely open a gap in the disc. The following gap opening condition is used (Crida et al. 2006): where RH = a · q 3 is the satellite's Hill-radius, cs the local speed of sound, q = M sat M planet and Re is the Reynolds number. As the gap opening in particular is very dependent on the choice of the scale-height h = cs Ω k , we can justify this prescription by recognizing that the disc is not selfgravitating (see Sect. 2.2.1) and we can assume that most of the disc mass will be concentrated around the midplane. However, recent works (Malik et al. 2015;Mueller et al. 2018) have shown that for massive discs, such as the one assumed here, the above criterion leads to easier gap opening compared to the results of hydrodynamical simulations. Moreover, Mueller et al. (2018) showed that one can better match the gap opening conditions in simulations by additionally requiring that the satellite's crossing time, i.e. the time it would take the satellite to migrate across a region with the width of the gap, τcross = RHS( da dt ) −1 (Malik et al. 2015) with RHS = 2.5 · RH (Pardekooper & Papaloizou 2009), is longer than the viscous time of the gas (τvisc = a 2 ν ), so Type II migration happens if: P < 1 and τvisc < τcross Once a gap is opened, the satellites migrate with the gap on the viscous time-scale given by: The factor bII =

1+
M sat 4π·a·Σgas (Syer & Clarke 1995) reflects that the migration slows down the heavier the satellite is compared to the gap.

Resonance trapping
Two satellites can get captured into resonances with each other, which tend to make them migrate together. This can be modeled such that as satellites move towards each other, there is a repelling force between the two. This will lead to an increase in the satellites' separation, which is modeled in our work like (Ida & Lin 2010): with RH = , where M planet refers to the planet around which the satellites orbit, is the Hill-sphere of their combined system. The inner satellite will feel half of this separation increase towards the central planet, while the outer one will feel the other half outwards.

Collisions
One of the restrictions of 1D models is that the collisions are trivial. To avoid an over-estimation of the satellite collisions (and therefore the final moon masses), we use a 2D approach when we test for collisions. Once a satellite is created, a random angle between [0, 2π] is assigned. This is updated with every timestep using ΩK and the test for collisions uses the actual 2D distance between two satellites. Hence, a collision happens only if this distance is smaller than the radius of the combined Hill-sphere of the two moons. When a collision happens in our model, we assume only the more massive one survives, and we add the mass of smaller body to the larger one.

Accretion
The satellitesimals accrete mass while migrating. The solid-accretion model by Greenberg et al. (1991) is used: where Rsat is the satellite's radius. The Σ dust in this case is the average dust density over the entire feeding zone of the moons. Greenberg et al. (1991) gave a value for the radius of this zone of R feeding = 2.3 · RH. This accretion is again multiplied by the same factor for the partial gap as the migration, because the dust and gas are assumed to be well coupled.
2.2.2.6 Dust depletion and refilling As the satellites are growing, the accreted mass is subtracted from the disc dust density, so the mass Mi in a given cell i decreases in the following manner: Like in the case of core-accretion simulations, the CPD is continuously fed from the circumstellar disc (Szulágyi et al. 2014), even in these GI simulations. We calculate the net mass influx from the SPH simulation:Ṁin,0 = 7.44 · 10 −5 M Jupiter year at t=0. Since the circumstellar disc will eventually dissipate, the feeding will decrease in the same way. Therefore, the feeding rate is also decaying exponentially on the disc dispersion time-scale:Ṁin =Ṁin,0e This influx has the same dust-to-gas ratio as the disc, in order to be consistent. As this is the total influx into the disc, a model is needed to determine how the disc reacts to changes in the dust density, specifically how it refills the dust gaps left by the accreting satellites. For that, a model like Cilibrasi et al. (2018) is used. It is assumed that the disc is in equilibrium at the beginning and that it wants to return to that equilibrium on a certain time-scale t refilling according to the following formula: This refilling can only ever be as much as the influẋ Min.

Varying Initial Parameters
Some of the initial conditions are a priori unknown, such as the dust-to-gas ratio. The following parameters are varied randomly in the different population synthesis runs: (i) The dust-to-gas ratio of the disc is varied in the range [10 −3 , 10 −1 ], based on observational constraints of circumstellar discs (Ansdell et al. 2016).
(ii) The dispersion time-scale of the disc is in the range [10 4 , 10 5 ]. The time-scales are intentionally short because gravitationally unstable CSDs undergo strong viscous transport as a result of gravitoturbulence (Durisen 2007), hence they are expected to dissipate faster than the typical Myr time-scale inferred for low mass T Tauri disc samples (see e.g. Mamajek (2009)).
(iii) The refilling time-scale of the disc is in the range [10 2 , 10 5 ]. This is the widest range possible (where it is still always bigger than the timestep), as this process is the least known.
(iv) The positions in which the initial embryos are created are randomized between 1% and 80% of the disc radius.
(v) The amount of initial satellite seeds are set between 5 and 20. It was tested that adding more seeds is not influencing the results.
The dust-to-gas ratios, refilling time-scales and initial satellite positions are distributed according to a log-uniform distribution, while the dispersion time-scales are distributed exponentially (Fedele et al. 2010).
Instead of randomizing all the parameters, one can also choose to set a fix value for a certain parameter. By doing this for different values of the fixed parameter, it is possible to investigate the parameter's influence on the results. This is done for the dust-to-gas ratio, and for the dispersion-and refilling time-scales.

CPDs at different distances from the star
We also investigate the effects of the planet's semi-major axis, i.e. the CPD distances from the star. Other than the 50 AU nominal case, we consider 5, 20, 35 and 70 AU semimajor axes. For this purpose, we rescale the CPD properties appropriately: (i) Because the CPD size depends on the planet's Hillradius (Quillen & Trilling 1998) -which linearly scales with the semi-major axis -the CPDs are scaled accordingly with 50AU a , where a is the planet's distance from the star. (ii) It is assumed that the CPD density profile is only dependent on the process that creates the planet (in this case gravitational instability) and the size of the disc (where the semi-major axis of the planet is relevant). The surface density profile is also assumed to be roughly the same regardless of the semi-major axis of the planet. This means that the 50 AU disc is mapped in such a way that the inner and outer edge match and the in between values are distributed relatively to the radii.
(iii) For the CPD temperature, a simple scaling is assumed, based on the distance of the planet from the star: a −1 .
(iv) The mass influx is also changed. Results from hydrodynamical simulations in Szulágyi (2017) showed that the influx decreases by a factor of 2 if the planet is moved 10 times farther out (so in this case the influx at 5 AU is double that of the influx at 50 AU), so the influx values are scaled according toṀin(a, t) = 3.25 · ( a AU ) −0.3 ·Ṁin(50 AU, t), where a is the semi-major axis of the planet.

Observational predictions
We also want to make a prediction on the observability of the satellites that are created in our population synthesis. For this, exoplanets that have already been discovered are considered, the satellites created in the population synthesis are placed in orbit around them (adjusted relative to the exoplanet size), and a detection probability is computed for every satellite, and averaged over all satellite systems. This way, it is possible to get an average probability per planet of detecting at least one exomoon. First, the detection probability for a given satellite around an exoplanet that orbits a star has to be determined. It has to be checked whether a given satellite is big enough compared to the star for modern instruments to be able to detect them via transits. The measure for this is called the transit depth, which is given as R sat R star 2 . If this depth is bigger than the instrument threshold, we set a detection probability of 1, 0 else. Second, geometrical factors need to be accounted for. If the inclination of the exoplanet's orbit relative to the observer is high enough, there is a small portion of time where an orbiting satellite will pass above/below the star and we will not be able to detect it. The fraction of time a satellite spends in a region where we can not detect it (and thus the probability that we can detect it) can be calculated with the following formula: where i is the exoplanet orbit inclination. This factor is used if cos(i) >= R star a planet +a sat , otherwise it is simply 1. Another geometrical factor is the time a satellite will spend directly in front of or behind the exoplanet and this factor we calculate like: and it is be used if cos(i) <= R planet a sat , else it is 1. These factors will then give a probability pj of detecting a satellite j around a given exoplanet. We then compute this probability for every satellite in a given satellite system created by the population synthesis and calculate the probability of detecting at least one of these satellites around a given exoplanet: This is done for every satellite system for a given exoplanet and averaged over the total amount of satellite systems to get the probability of detecting at least one satellite around a given exoplanet and then for every exoplanet in the catalogue (where the necessary information is found) to get in the end an average probability of detecting at least 1 satellite around an exoplanet. Since the exoplanets considered have varying masses and orbital radii, the satellite properties are rescaled to better fit around those planets.
(i) Depending on the exoplanet's orbital radius, the satellite systems from population synthesis run with the CPD at a specific semi-major axis are considered. So if the exoplanet has an orbital radius < 13 AU, the satellites from the 5 AU case are considered, if the exoplanet is in the range 13 AU < a exoplanet < 28 AU the results from the 20 AU case are considered etc.
(ii) The satellites' orbital radii are then rescaled as follows: where R CPD planet is the radius of the planet around which the CPD was formed and B is a constant defined such that a sat, rescaled = R exoplanet if a sat, original = R CPD planet .
(iii) The satellite masses are assumed to be directly dependent on planet mass, so that the masses can be rescaled like where m CPD planet = 10 MJupiter is the mass of the planet around which the satellites from the population synthesis are created.

RESULTS
We present our results in the following way. First we analyse a reference case, that of a CPD at 50 AU from the central star, describing the properties of the emerging satellites' population, and how those vary as we vary the key parameters such as dust-to-gas ratio, disc dispersion time-scale and dust refilling time-scale. Second, we compare the same properties for CPDs forming around disc instability proto-planets located at different distances from the central star.
3.1 The reference case: a CPD at 50 AU satellite is created, with a maximum of 10 in less then 1 % of the simulations. In the 35 % of cases that no embryo reaches the cut-off mass, there are 0.045 % cases where there is not even an embryo alive at the end of the simulation, which means they all either migrated into the planet or are lost to collisions. The dependence of the number of surviving satellites on dust-to-gas ratio and refilling time-scale is as expected. It reflects the fact that it is easier to reach a mass higher than 10 −5 M planet when there is more dust, or when the gap is refilled faster. Varying values for fixed dispersion timescales give a non trivial picture though. Fig. 5 shows that as the dispersion time increases, the number of satellites first drops, then increases again. The reason for this behaviour is that the first generation of embryos will be the one that accretes the most mass, as the disc will have more dust available at earlier times. If the dispersion time-scale (and thus the disc lifetime) is low, less of those first generation embryos will eventually migrate into the planet. For intermediate dispersion time-scales, many first generation satellites will be lost into the planet, but subsequent generations will not have the time to grow enough to replace the population of high mass satellites. This can only happen in long dispersion time-scales, which is the behaviour Fig. 5 shows. Fig. 6 shows the distribution of satellite masses as a fraction of the mass of the central planet, namely in units of 10 Jupiter masses. We see that many are in the range of Jovian satellites, between 10 −6 and 10 −5 planet masses, but they can reach as high as 10 −3 planet masses, which corresponds to 3 Earth masses. Interestingly, such natural upper mass limit emerging from our model is of the order of the mass of observed Super-Earths, at and above which planets with H/He envelops are known to exist, suggesting that following the evolution of the upper end of our mass function might Figure 6. The distribution of the masses of all satellites that survive until the end of the simulation. The masses below 10 −6 M planet are cropped because of the way we spawn embryos, there are a lot of them that get created late enough that they simply do not accrete mass anymore, so the numbers are heavily skewed towards them. require to include the possibility that they could accrete a gaseous envelope. Varying values for dust-to-gas ratio and refilling time-scale yield a trivial outcome; higher dust-to-gas ratio implies higher masses, and shorter refilling time-scale implies more satellites can form at any given mass. The results for the dispersion time-scales are shown in Fig. 7. The behaviour corresponds to the result in Section 3.1.1; for fast dispersion times, there are many high mass satellites, close to 10 −3 planet masses. These satellites vanish as the dispersion time-scale becomes higher, because they migrate into the planet. As the dispersion time-scale increases, there are more satellites in general, with a similar distribution, as now the subsequent generations have more time to accrete mass, but do so on a slow enough time-scale that they hardly migrate into the planet in the remaining available time. Fig. 8 shows the integrated mass of the different satellite systems created. In our Solar System's gas giants, the satellite systems are consistently around 2 · 10 −4 planet masses. Figure 8. The integrated mass of the satellite systems. There is a peak around 10 −5 planet masses, which is a little lower than what we observe in our Solar System's gas giants. It also shows a wide spread, with the majority between 10 −5 and 10 −3 M planet Figure 9. This plot shows the amount of mass lost into the central planet against the combined mass of surviving satellites. It shows that retained mass depends almost linearly on lost mass and that the retained mass is lower or about the same as the lost mass About 50 % of the cases are between 10 −5 and 10 −4 planet masses, lower than what can be observed in the Solar System. About 30 % are 10 −4 planet masses and higher. As the maximum integrated mass is not much higher than the maximum mass possible for a single satellite, it also shows that satellites that massive are very rare and tend to dominate their systems.

Mass distribution
Previous work (Bolton et al. 2017) has shown that Jupiter is about 20 times more enriched in metals relative to the solar metallicity value. Pollution by embryos migrating into the planet could explain this result. Fig. 9 displays the mass retained by the satellites versus the mass lost by them, whereby the median mass lost is about 3 · 10 −3 planet masses, or about 10 Earth masses. It shows that mass that is lost scales with the retained mass, where the latter is always lower than the former.
The results for varying values of fixed parameters mirrors those in Section 3.1.1 and of the mass distribution: higher dust-to-gas ratio and faster refilling time-scales results in higher lost and retained masses. Fig. 10 shows that as dispersion time-scale increases so does the lost mass, which corresponds to the decrease in the number of heavy satellites which get lost into the planet. This also accounts for the decrease in the maximum retained mass. For long dispersion  Figure 11. Positions of the satellites with at least Europa mass at the end of the simulation. It shows a peak near the planet, with orbits a small multiple of the planet radius, which is what we would expect. However, it also shows another peak in the outer disc, which is more unexpected and a result of low gas density, which makes satellites stay in the region they are created. This allows them to accrete mass throughout the simulation without migrating into the planet. time-scales, both lost and retained mass increases, as then subsequent generations of satellites can accrete more mass. Fig. 11 shows the distribution of final positions for satellites with at least Europa mass. There is a peak close to the planet, where the semi-major axes of the orbits are comparable to the planet radius. Even though we account for resonant trapping, there are no resonant configurations that show up consistently, as they would manifest as separate peaks. As the distance from the planet increases, there are less satellites orbiting. Migration time-scales scale mainly with gas density and satellite mass and as the feeding zone will be small this close to the planet, satellites will have to accrete mass further out and migrate into this region. As gas density increases towards the planet, their migration will Figure 12. Formation temperature of the surviving satellites of at least Europa mass, which shows that most of them will form at temperatures below the freezing point of water and thus should all be icy satellites.

Radial distribution of satellites
speed up, even though their mass will remain very stagnant and thus it is much more likely that any satellite entering the region past the close peak will migrate close to or into the planet eventually, leaving an evacuated area.
It is more surprising that a second peak is located very far out. The expectation would be that a satellite which reaches that much mass would migrate away from its formation site so that the distribution should spread out more than what the plot shows. The behaviour of satellites that are far out is investigated by specifically placing embryos in the region past the far peak and monitoring their evolution. These tests show that only embryos that accrete mass very fast in their early history can travel past the far peak. This is because migration speed is proportional to both gas density and the mass of the migrating object. The gas density (which in the region past the far peak is already very low) will drop over time because of the dispersion of the disc, but the mass will increase as the embryo accretes mass. For a satellite to be able to migrate past the far peak, it needs to overcome the low gas density and its drop over time with a fast enough accretion early in their life. This can only happen with certain combinations of dust-to-gas ratio, dispersionand refilling time-scales. For most satellites, their accretion is too low to overcome the drop in gas density, resulting in very little migration. This means however that they can accrete mass throughout the whole disc lifetime, undergoing very few collisions with other proto-satellites, which makes it possible for them to reach the mass of Europa or higher. Therefore, the peak on the outer edge is a result of the interplay between satellite accretion and gas density.
We also analysed the changes in the results as we vary the main initial conditions separately. For higher dust-togas ratio, more satellites migrate close to the planet, because they grow more massive, which facilitates migration. The refilling time-scale has the same effect, although less pronounced. For short dispersion time-scales, there are also more satellites closer to the planet. These are the high mass satellites we discussed earlier, and as the dispersion timescale increases, these will migrate into the planet, resulting in a shrinking peak in the histogram.

Formation temperature
We also look at the temperature of the gas at the starting location of each satellite, as that should be indicative of the composition of the satellites. As water in the disc freezes at about 180 K (Lodders 2003), we can at least use the formation temperature to estimate the likelihood to obtain icy or rocky satellites. Fig. 12 shows that the vast majority of satellites form in regions with temperatures below 100 K, which is a strong indication that in discs like this, most satellites should end up being icy. This is not surprising, as Fig. 3 shows that, even at the beginning, most of the disc is near or below 180 K.
We again checked the influence of varying the initial parameters one by one. The refilling time-scale has no apparent influence on the shape of the distribution. Fig. 13 shows the result for three different values for the dispersion time-scale. We see that for short dispersion time-scales, the formation temperatures are very high. That is because the short disc lifetime, which is coupled to the dispersion time, means that the surviving satellites will be older, as it is hard for them to migrate into the planet, and thus are created in an earlier, hotter disc. In Fig. 14 we investigated the effect of the dust-to-gas ratio. It shows that for low dust-to-gas ratios, the satellites will also be created in a hotter disc. That is because for this low amount of dust, the only embryos that will grow significantly in size, are necessarily older.

Time-scale to reach the Europa mass scale
We also asses how fast the satellites form and use the time in which they reach Europa mass as a scale. In principle, their final composition and structure can be heavily influenced by how fast they form, and subsequently how much time they spend in the gas disc as satellites. We choose the mass of Europa as a reference because that corresponds 10 −6 the mass of the parent planet in our case, which turns out to be the lowest mass significant satellites will have. The  majority of satellites will reach the Europa mass scale on a time-scale similar to the dispersion time-scale, between 10 4 and 10 5 years, which is quite long. The maximum is around 10 6 years, which is just shorter than the maximum time a simulation can run in general (14 · 10 5 years). However, it can also be very short, 10 3 years and slower, although such cases are rare.
We also consider the influence of the other other free parameters on the latter time-scale. Variations of the dustto-gas ratio produce the expected result: the higher the ratio, the more satellites reach Europa mass very fast, which results in a much wider distribution. If we fix the dispersion time-scale, the distribution looks essentially the same as in the fully randomized case, but, since the disc lifetime is coupled to the dispersion time-scale, the bulk of the distribution shifts along with the dispersion time-scale. In Fig. 16 we fixed the refilling time-scale at various values. It shows that the refilling time-scale only influences the very fast growing satellites, as a fast refilling means more satellites with a Europa mass time-scale in the 10 3 years region. The rest of the Figure 16. Comparison of Europa mass times in systems with a set refilling time-scale, but with the other parameters still randomized. Figure 17. This plot shows a comparison of how many satellites with a mass greater than that 10 −5 M planet are surviving for the different radial distances of the CPD from the star. distribution is largely the same, which means that the satellites with long Europa time-scales are satellites that migrate consistently enough, so that their mass comes from the disc itself, and not through influx into the gaps. Fig. 17 shows a comparison of how many satellites with masses higher than 10 −5 planet masses survive until the end. Figure 18. This plot shows a comparison of the masses of surviving satellites for the different central planet positions. We see that as the planet is further out (and thus has a bigger disc), the satellites become in general heavier.

Number of satellites
The results are very similar, there is only a very slight increase in numbers as the disc grows in size. This is the result of two opposing effects: a bigger disc (which means a bigger semi-major axis of the planet) has more space for the satellites to transverse and more mass for satellites to accrete, so with increasing size the satellites should be more massive and thus more numerous. But as the distance from the star increases, the influx into the disc decreases, which means that embryos on average grow slower, which makes satellites less numerous. These two effects in this case result in the observed behaviour. We find that for the 5 AU and 70 AU cases, there is always at least one embryo that survives until the the end of the simulations, even if it doesn't reach 10 −5 planet masses. For the 20 AU, 35 AU and 50 AU we find that there are 0.045 %, 0.045% and 0.02 % cases respectively where no embryo survives until the end. Fig. 18 shows the mass distributions. It shows that as the disc grows in size (and thus also in mass), the average mass of a satellite increases as well, as shown in the maximum masses reached. In the 5 AU case, satellites grow at most to about 1-2 Earth masses, while in the 70 AU case the maximum is at about 10 Earth masses. Secondly, it also shows that lower mass satellites grow relatively less likely, while higher mass satellites become more likely. In the 70 AU case, this even results in a slight peak close to the maximum, which is not only because there is more mass in the disc, but also that satellites are less likely to be lost into the planet and so more first generation satellites should survive. The integrated mass shows corresponding behaviour: Figure 19. The integrated mass of the satellite systems for the different radial distances of the planet. They all have a maximum between 10 −5 and 10 −4 planet masses, but as the disc grows bigger the distribution grows wider.

Mass distribution
the increasing average mass means the distribution grows wider towards high masses. Meanwhile the peak moves towards lower mass and shrinks, as some of those satellites grow heavier in a bigger disc.

Radial distribution of satellites
In general, the 20 AU to 70 AU cases show very similar behaviour: a peak very close and a peak very far away and an evacuated middle section. The percentage of satellites close to the planet increases as disc size increases, because, as they on average migrate longer, it is easier for them to stop close to the planet instead of migrating into it. However, the percentage of satellites far out also increases, because with a bigger disc, the outer parts have more mass, so it becomes easier for an embryo to accrete a lot of mass. The middle section at 20 AU also shows a slight peak which vanishes as the disc size increases. As the disc grows in size, it will be harder for a satellite to move past the far peak and into this middle region. This is not only because the effective distance to migrate increases, but also because of what we established in Section 3.1.3, where we found that a satellite will need to accrete a lot of mass very early in its life in order to be able to migrate past the far peak. This becomes harder as the disc grows because the mass influx decreases and since satellites this far out tend to stay in a small region, because of the low gas density and migration, this means that their accretion time-scale is closely related to how fast the mass in their feeding zone is replenished. The 5 AU case is somewhat different than the other cases. The histogram of satellite locations does have a peak far out in the CPD, like the other cases, but it does not show a peak close to the planet. Instead, the peak in the middle region is Figure 20. This plot shows the amount of mass lost into the central planet against the combined mass of surviving satellites. It shows that the retained mass depends almost linearly on lost mass and that usually the mass retained is lower or about the same as the lost mass. Figure 21. Comparison of satellite locations at the end of the simulations for the varied semi-major axis cases. When the planet is at 20 AU or further out, the histograms look roughly the same as the nominal (50 AU) case: one peak of satellites in the inner CPD and one peak in the outer disc. As the disc grows bigger, the inner peak grows as well, showing that as the disc size increases, it becomes less likely for satellites to migrate into the planet. Figure 22. Comparison of formation temperature of the surviving satellites for the different cases. We see that even though the disc becomes hotter the closer the central planet is, most satellites will still be formed icy. For the 5 AU case, the peak at 1500 K is because at this temperature the dust evaporates, so there are no embryos formed at temperatures above. more prominent. This is an effect of both the small disc size and the higher mass influx. The influx makes it easier for embryos to accrete mass fast and then migrate past the far peak. However, as the disc is small, embryos and satellites alike hardly stop when they are in the vicinity of the planet, rather will migrate into the planet more easily.

Formation temperature
The formation temperature shows no big surprises. As you get closer to the host star, the disc heats up and thus the average temperature at which the surviving satellites are formed increases as well, albeit only slightly for the majority. It still shows that until 20 AU, most will form icy. The 5 AU case is a bit different, as this disc becomes hot enough in certain sections that dust evaporates. Thus there is a peak in formation temperature at around 1500 K, where dust evaporates. These will form late enough into disc lifetime that migrating into the planet will be very difficult. It also shows that in this case the satellites form significantly hotter in general.

Time-scale to reach the Europa-mass
For all the cases, most of the satellites that reach the Europa mass scale do so between 10 4 and 10 5 years. However, the distribution widens towards fast time-scales as the disc grows bigger. There are two effects that play into this: first, since bigger discs have more mass, the average mass in the feeding zone of an embryo is higher, which allows for faster Figure 23. Comparison of time-scales to reach a mass equivalent to that of Europa in the different cases. Most of them peak around 10 4 years, but as the disc grows bigger, the distribution widens. This is because in general the cells which make up the disc will contain more mass and it will thus be easier for satellites to accrete mass very fast. accretion. Second, as you get closer to the star, the smaller disc size means embryos that accrete mass very fast are more likely to migrate into the planet. We also see that the peak of the distribution moves towards longer time-scales as the disc grows bigger. That is because on average, the embryos will form farther out in the disc, so their orbital time-scale increases and thus their accretion time-scale increases as well. This means that, on average, the farther out an embryo is, the longer it would take it to reach Europa mass. Fig. 24 shows a plot of transit depths for every satellite created in the population synthesis with a mass greater than Europa. If they are smaller than that, they have a transit depth so low that no instrument could detect them, so we left those out of the plot. Even then, most of the satellites are too small to be detected for any instrument but E-ELT (before taking into account the geometrical factors). This already shows the difficulty of detecting exomoons such as those generated in the population synthesis. The averaged probability to find at least one satellite around an exoplanet confirms that, as shown in Table 1. Even with the best instrument, the probability of find- Figure 24. Plot of the transit depths for the satellites with at least Europa mass (any smaller ones will have a transit depth below all the equipment thresholds). The colored lines represent different instrument thresholds, therefore the population lying to the right of each line is that could be theoretically observed by the given instrument.

Observational Predictions
ing an exomoon around a given exoplanet is only around 3%. This is not only because the satellites are of small mass, and hence size, but also because the heavier satellites tend to orbit close to the planet, which means that they spend more time in front of or behind it. For the satellites farther out it is the opposite; they tend to be smaller, but they also spend less time obscured by the planet.

DISCUSSION
In the population synthesis approach several assumptions are needed to design a modeling framework simple enough to be flexible and allow studying the role and interplay of the various underlying physical processes. Such assumptions have a varying degree of impact on the results, which we discuss here. The first choice we made concerned the evolution of the disc. We started from gas density and temperature profiles that resulted from hydrodynamical simulations of Szulágyi et al. (2016b) that generated a 10 MJupiter protoplanet and its CPD at 50 AU. We employed a simple model to evolve both density and temperature of the CPD exponentially on a certain time-scale, neglecting completely local thermodynamical effects. We also reduced the 3D disc structure extracted from the hydro simulation to construct a 1D model of the disc by using azimuthally averaged quantities, in accordance with Pringle (1981). Parameters such as scale-height, surface density and sound speed, extensively used in our model, are all born out of this simple 1D approximation.
Another choice we made concerns the generation of embryos. We used a seeding approach, where we set a certain set of randomly distributed locations at which we inserted the embryos. As a result, the creation of embryos was not tied to the disc parameters, apart from, implicitly, to the disc size.
As we investigated the influence of the planet's semimajor axis, we rescaled the CPD properties (temperature, density) accordingly. For the temperature, we prescribed that the CPD has a higher or lower temperature normal-ization depending on the distance of the parent planet from the star. We also assumed that the form of the surface density profile does not change irrespective of the semi-major axis. We know that disc instability tends to generate gas giants from fragmentation predominantly in outer region of the CSD (at R > 30 AU), so if a planet is closer to the star, the more likely it e.g. migrated there from the CSD outer regions. The migration of the CPD, and the variation in temperature due to the varying stellar irradiation flux at varying distances would affect the density-and temperature profiles of the CPD significantly.
There are also some parameters that we do not vary for simplicity, while we expect they should vary. One obvious example is the initial seed mass, which we set at 10 −8 planet masses. We tested varying that parameter, and determined that even lowering to as much as 10 −10 planet masses has no significant influence on the results, because in the first few accretion cycles the embryos grow quickly to 10 −8 M planet .
We also set the seed locations at the beginning of the simulation and created later generations at the same locations instead of then varying it. We tested that using a different approach, namely randomly choosing the spawning locations at each subsequent iteration, does not affect the results.
We also set the number of seeds between 5 and 20. We tested that decreasing or increasing this number does not have a big influence on the results, because it is extremely rare to have more than 10 massive satellites in the CPD at a given time.
It is noteworthy that there is no other work on satellite formation in disc instability scenarios, nor works that address in general satellite formation around planets much more massive than Jupiter as in the present case, hence we could not compare directly with any previous work. However, we can still attempt a qualitative comparison with previous literature, which has exclusively focused on satellite formation in core-accretion and around a Jupiter sized planet, to address the formation of Galilean moons. Miguel & Ida (2016) used a minimum mass disc based on how the Jovian satellite system looks today with a cavity between disc and planet. A seeding approach was also used and accretion and migration of the embryos were taken into account, as well as dust and gas evolution. Temperature was not evolved and no dust influx from the stellar nebula was considered. The results showed a similar mass distribution as in our work, although with lower masses in general due to the disc being smaller and lighter. It also seemed to mirror our results for the positions of the orbits, with a cluster close to and a cluster far from the planet, although this may be influenced by the gas density which had a steep incline in the middle regions. The cavity also might influence how many satellites end up very close to the planet. Cilibrasi et al. (2018) used a model very similar to ours in terms of conceptual framework, but they developed that in context of a CPD formed by core-accretion, also drawn from a 3D radiative simulation as in our case. Apart from the different CPD structure emerging in different formation scenarios, in the latter work a dust evolution model was used to determine the initial dust density, which naturally produced a region with increased density, a so-called dust trap. The satellites were assumed to form there through streaming instability, which requires the local dust-to-gas ratio to be greater than unity. This alone leads to a very different mass distribution, very sharply peaked unlike the wide distribution in our work. That is because the satellite formation was tied to disc parameters. As a result it had to stop at a certain point, and since the satellites accreted most of their mass from the dust trap, very few remain close to the initial mass. Fujii et al. (2017) used a more restrictive population synthesis to take a detailed look at the orbital evolution of moons. They evolved both surface density and temperature of their disc (around a 0.4 MJupiter planet) by numerically solving appropriate diffusion equations. They allowed their satellites to migrate and use a more extensive model of resonance trapping, and their satellites started with a fixed mass comparable to Io. They varied both influx and viscosity to test at different scenarios. Their results showed that, with the inclusion of proper resonance trapping, they could almost always get out a 1:2:4 resonance, that is present in the Galilean system. The resulting satellite positions ended up different in their different runs, but the resonance was very consistently present. This is opposed to our model where we found no consistent resonant configurations, probably because of the simpler treatment of resonance capturing.

CONCLUSIONS
In this work we investigated the formation of satellite systems around a 10 MJupiter gas giant created by gravitational instability with a semi-major axis of 50 AU. We used the gas density and temperature profiles as well as the mass influx of a disc obtained with an SPH simulation (Szulágyi et al. 2016b) and assumed that the dust had the same profile as the gas, but multiplied by the overall dust-to-gas ratio of the disc.
We then used a population synthesis approach, where we seeded the disc with embryos. The CPD itself was evolved over time on the dispersion time-scale and refilled with dust on a refilling time-scale. The embryos migrated, accreted mass and collided while also allowing for subsequent generations to be formed.
We randomized some of the not well constrained initial conditions, such as the dust-to-gas ratio, the disc dispersionand refilling time-scales as well as the number and initial positions of the seeds at the beginning of every run. We also investigated the influence of the semi-major axis of the CPD on the satellite systems by rescaling the size, temperature and mass influx of the disc accordingly.
In the nominal case (planet with a CPD at 50 AU from the star), our results showed that in ∼ 60 % cases, 1-4 satellites formed with a mass of minimum 10 −5 planet masses, similarly to what we can observe in our Solar System's satellite systems. The mass distribution peaked at around 10 −3 planet masses (i.e. 3 Earth masses), with most of the satellites having had masses in the range of the Jovian satellites. The integrated mass had a peak at around 2 · 10 −5 planet masses, an order of magnitude lower than what we can observe in the case of Jovian, Saturnian and Uranian moon systems. Migration of the satellites led to ∼10 Earth masses of solids lost into the planet. We also found that the vast majority of satellites were created at temperatures much lower than the freezing point of water and as such are very likely to be icy. This is no surprise, as the temperature of the disc was very low to begin with. Most satellites also formed between 10 4 and 10 5 years, comparable to the dispersion time-scales, with some having formed very rapidly in around 10 3 years.
The influence of the semi-major axis of the planet on the number of satellites and their formation time-scales turned out to be minimal. The mass histogram of moons showed that it became harder for heavy satellites to form in CPDs closer to the star. First, because they could easily migrate into the planet and second, because the overall CPD mass was smaller to make moons.
We also found that the probability of detecting satellites like the ones in this population synthesis is very low (≤ 3%) even with E-ELT. This is mainly because the satellites are expected to form beyond the snowline (around giant and ice giant planets), and the current detection methods are not sensitive to this further-out population. On top of that, the heavier satellites tend to orbit close to the planet which makes their detection more challenging. However, giant planets can get closer to the star via various dynamical evolution processes, and if they could keep their moons during this process, then these satellites should be easier to detect.