An N-body population synthesis framework for the formation of moons around Jupiter-like planets

The moons of giant planets are believed to form in situ in Circumplanetary Discs (CPDs). Here we present an N-body population synthesis framework for satellite formation around a Jupiter-like planet, in which the dust-to-gas ratio, the accretion rate of solids from the Protoplanetary Disc, the number, and the initial positions of protosatellites were randomly chosen from realistic distributions. The disc properties were from 3D radiative simulations sampled in 1D and 2D grids and evolved semi-analytically with time. The N-body satellitesimals accreted mass from the dust component of the disc, interacted gravitationally with each other, experienced close-encounters, both scattering and colliding. With this improved modeling, we found that about $15\%$ of the resulting population is more massive than the Galilean one, and only $8.5\%$ were in resonances. The moons reach Europa's mass most frequently in $10^5$ years. In $10\%$ of the cases, moons are engulfed by the planet, and $1\%$ of the satellite-systems lose at least 1 Earth-mass into the planet, contributing to the giant planet's envelope's heavy element content. In $1\%$ of cases, we found eccentricities and inclinations of moons above 0.1. We examined the differences in outcome between the 1D and 2D disc models and used machine learning (Randomized Dependence Coefficient together with t-SNE) to compare our population with the Galilean system. Detecting our population around known transiting Jupiter-like planets via transits and TTVs would be challenging, but $14\%$ of the moons could be with an instrumental transit sensitivity of $10^{-5}$.


INTRODUCTION
Two main theories of giant planet formation in Protoplanetary Discs (PPDs) have been developed in the last few decades, the Core Accretion Model (CA; Pollack et al. 1996) and the Gravitational Instability Model (GI; Boss 1997). The Core Accretion scenario assumes that dust grains orbiting in a PPD are able to coagulate and grow efficiently into boulder-sized objects above meter sizes, thus overcoming fast drag-induced inward drift towards the star (Rice et al. 2004;Lambrechts & Johansen 2012;Drążkowska et al. 2013), and then later into kilometer-sized planetesimals, which can then form planetary-sized cores by gravitational accumulation (Lissauer 1993;Levison et al. 2015). When the cores are massive enough, they start accreting gas, until they achieve the final structure of a giant planet composed of a solid nucleus and a gaseous envelope. In the Gravitational Instability (GI) scenario (Boss 1997;Durisen et al. 2007), giant planets originate from gas clumps that form in a gravitationally unstable PPD by direct collapse. These possibly accrete solid material only in a later phase, eventually potentially developing solid cores and enriching their envelopes (Helled et al. 2006;Helled & Schubert 2008;Boley et al. 2010).
In the final phase of gas giant planet formation, a Circumplanetary Disc (CPD) made of gas and dust develops around the nascent planet (e.g. Kley 1999;Lubow et al. 1999;Canup & Ward 2002;Ayliffe & gravitational potential well (compact solid core vs. collapsing gas cloud) .
Efforts have been made also in observational astronomy, with both attempts to directly detect of CPDs or even exomoons and predictions for observability stemming from hydrodynamical numerical simulations (Szulágyi et al. 2018a;Isella & Turner 2018;Szulágyi et al. 2019;Szulágyi & Ercolano 2020) and from modeling features of moon transits such as transit timing variations (TTVs), transit duration variations (TDVs), and apparent planetary transit radius variations (TRVs) (Kipping 2009;Rodenbeck et al. 2020). A first observational evidence of the presence of a CPD was found by multiple groups in 2019. Christiaens et al. (2019) found infrared excess due to circumplanetary material with VLT/SINFONI around the ∼ 10 Jupiter-mass protoplanet PDS70b. Isella et al. (2019) detected continuum emission with ALMA near two protoplanet candidates. They were able to identifying one of them as a potential CPD with 2 − 4.2 × 10 −3 ⊕ dust-mass (assumed total mass of 10 −5 − 10 −4 planetary masses). Regarding satellite system detection, after the first indirect evidence of V1400 Centauri (Kenworthy & Mamajek 2015), where gaps in a detected ring system were supposed to be caused by exomoons, the first candidate has been identified as a Neptune-size moon orbiting the 10 Jupiter masses planet Kepler-1625b(Teachey & Kipping 2018. Circumplanetary Discs are considered the birth places for satellite systems around giant planets (Canup & Ward 2002). As we see in our Solar System, giant planets often have systems of regular satellites, massive enough to be spherical, with low eccentricities and inclinations, similar to a scaled-down version of a planetary system. Due to the so-far gas only hydrodynamic simulations of the circumplanetary discs, we do not know yet what dust content they might have, or what grain-size distribution one can expect. It is known that the ratio between the total mass of moons around Jupiter, Saturn and Uranus and their planetary mass is about 2 × 10 −4 , implying a minimum total gas mass of the CPD of 2 × 10 −2 planetary masses (computed with the interstellar medium value of dust-to-gas ratio: 1%). This is called Minimum Mass Sub-nebula Model, or MMSM (Mosqueira & Estrada 2003). Such a massive disc would imply some difficulties for satellite formation, such as too fast migration rates (i.e. all the moons can be lost by migrating into the planet) or too hot temperatures for ice to condense, but a continuous feeding from the PPD through the meridional circulation (Szulágyi et al. 2014) would allow a much lighter CPD to exist over time, while still able to form enough satellites, in what is usually called gas-starved scenario (Canup & Ward 2006). The persistent flow of dust from the PPD could possibly allow the formation of more than one generation of satellites, migrating one after the other into the central planet, in the so-called sequential formation scenario (Canup & Ward 2002).
There have been several studies of satellite formation in Circumplanetary Discs so far. Fujii et al. (2017) used a 1D analytical model for moon formation and migration in a CPD, and found that moons are often locked in resonant chains, in agreement with the case of the Galilean satellites. Heller & Pudritz (2015) studied the evolution of a CPD also taking into account different heating processes, focusing on the position of the ice-line and concluding that satellites should form in the very last phase of a giant planet formation, when the environment is sufficiently cold for icy satellite formation. Other works have been focusing on the accretion processes of satellites. Other than the satellitesimal accretion and dust accretion, the pebble accretion scenario has been investigated in the last few years. Shibaike et al. (2019) showed that slow pebble accretion of planetesimals captured in CPDs would be able to prevent a too fast migration process and, given some fine-tuned halt conditions and physical constraints, to build Galilean-like systems in masses, positions (and resonances), and compositions. Ronnet & Johansen (2020) showed that pebbles would not be able to flow from the PPD to the CPD because of the pressure gradient at the gap edge, while the smaller grains are captured by gas drag. This mechanism would ablate pebbles, but providing the circumplanetary disc with small dust grains, which then coagulate into pebbles inside the CPD. This pebble accretion scenario furthermore represent rapidly migrating satellites, which are piling up at the inner boundary of the CPD, creating resonant chains.
A common approach in studying satellite formation is the so-called population synthesis method, that have been already widely applied to planet formation (Ida & Lin 2004;Benz et al. 2014;Mordasini 2018). Using such an approach means running many thousands of individual simulations, varying initial conditions, for example the moonlets' initial positions, the disc dust-to-gas ratio, its viscosity, or lifetime. An early example of an application of this method to moon-formation comes from Sasaki et al. (2010), in which the authors numerically computed the viscous evolution of a CPD, following the prescriptions of Pringle (1981), including a cavity between the disc and the planet, solid accretion onto the satellite-seeds and type I migration as in Tanaka et al. (2002). Their simulations resulted in systems of four-five satellites, often locked in resonances. A similar approach has been followed by Miguel & Ida (2016), in which the CPD was modeled as a MMSM, with randomly choosing the initial positions of 20 seeds, their initial masses, and the dust-to-gas ratio of the disc. Their results show that the high density gas causes very rapid migration of satellitesimals. Nevertheless, the conditions in the disc also favor the formation and the survival of large satellites, but massive and faraway bodies are difficult to reproduce. Another population synthesis approach was shown in Cilibrasi et al. (2018), where the CPD was an azimuthally and vertically integrated version of a 3D radiative hydrodynamics simulations of Jupiter's CPD (Szulágyi 2017). Up to 20'000 simulations were run with different dust-to-gas ratios, disc life-times and timescales of dust refilling from the PPD, showing a great occurrence of sequential formation and massive satellites in this 1D limit. Similar models have been also implemented to study satellite formation around ice giants in the Solar System (Szulágyi et al. 2018b) and around Jupiter-like exoplanets that are forming with Gravitational Instability scenario resulting in significantly different CPD characteristics (Inderbitzi et al. 2020).
Another approach to study satellite formation is using N-body simulations. Instead of solving the motion of satellites and their migration in a 1D scenario as all previous examples, a proper N-body integrator can reach a much better accuracy in the interactions between satellites (including resonances), as in Moraes et al. (2018). Those authors simulated a MMSM disc, in which embryos could move in a 3D space, following gravity interactions, together with smaller satellitesimals, that were the main source of accretion. They run a few different systems with the MERCURY N-body integrator (Chambers 1999), showing the influence of the number of initial seeds, the temperature of the disc and its density on the final architecture of the systems.
In this paper, we present an N-body population synthesis framework for satellite formation around Jupiter-like planets, combining the advantages of the last two described approaches. In our model, satellites were treated as N-body particles in an integrator including close encounters (Chambers 1999;Grimm & Stadel 2014) and individual time-steps (Saha & Tremaine 1994), while more than two thousand simulations were run with different initial conditions, i.e. with randomly choosing the dust-to-gas ratio, the solid accretion rate from the PPD and the initial positions and number of the seeds. We also included dust accretion onto the moonlets from a dust disc, and the interaction between the seeds and the disc itself, resulting in migration, eccentricity and inclination damping (Cresswell & Nelson 2008). Furthermore, following the approach of Cilibrasi et al. (2018), our model does not replicate a MMSM or gas-starved disc, but uses 3D thermo-hydrodynamical simulations from Szulágyi (2017) to build a 1D disc model, in which the N-body framework is embedded along with semi-analytical recipes to describe processes involved in the interaction between the satellitesimals and the surrounding CPD. We also tested the evolution of satellites, particularly their accretion, when the integrator is embedded in a 2D disc model, from Szulágyi (2017). We compared the results of the 1D disc + 3D N-body approach, and the 2D disc + 3D N-body case.

METHODS
The model presented in this paper is a semi-analytical model, consisting of a CPD in which satellites are treated as growing embryos that migrate through the disc and interact with each other. Since migration and accretion of protosatellites depend on the disc properties, the disc is also evolving, meaning that it is dissipating and cooling, while the dust profile is evolving accordingly to satellite accretion. In the model, the units are always for distances, for masses, for temperatures and for times.

Disc Structure
The circumplanetary disc has been modelled using the results of the simulations described by Szulágyi (2017) for Jupiter's circumplanetary disc. These simulations were the base also for our previous work (Cilibrasi et al. 2018), but this time using the CPD model with the coldest planetary temperature (1000 K instead of 2000 K as in Cilibrasi et al. 2018) among the ones presented in the Szulágyi (2017). As in Cilibrasi et al. (2018), the CPD has gas and dust density profiles, a temperature profile and the same parameters as in Szulágyi (2017): = 0.004 for viscosity (Shakura & Sunyaev 1973), the heat capacity ratio = 1.4 (molecular hydrogen), and the mean molecular weight = 2.37, that all have a role in the migration of satellite seeds. The planet is similarly assumed to be a Jupiter-equivalent at 5.2 AU, with = 1 and = 1 . In this model, the radial profiles of the CPD are built on a 1D radial grid, ranging logarithmically from 1 to 500 (about 2/3 of the Hill radius) with 1000 cells in total. A cavity between the planet's surface and the disc is not considered in the model, similarly to the original hydrodynamical simulations. This is because the magnetic field of the young Jupiter and the ionization of the CPD might have not been strong enough to produce such a cavity at the time of formation (Owen & Menou 2016). The initial profiles of the disc were taken as interpolations of the values coming from the hydrodynamic simulations, averaged both azimuthally and vertically in the case of the 1D density profile, or taken in the mid-plane as in the case of the temperature. Averaging azimuthally greatly decreases the computational expense and introduces a negligible approximation.
In fact, as explained in Section 2.3, the disc structure interacts with the forming satellites via two processes: the gas component via migration and the dust via accretion. The key parameters on which migration depends are the gas surface density Σ and the local aspect ratio of the disc = ℎ/ . These two quantities are defined as follows . Gas (blue) and dust (black) density profiles of the circumplanetary disc at the beginning of the simulations. The dust-to-gas ratio here is chosen to be 1% as an example, but varied between the different runs.
where is the gas volume density, = √︃ = √︃ is the sound speed in the mid-plane (with being the Boltzmann constant) and Ω is the Keplerian angular velocity. Furthermore, migration is always much slower than orbital motion (depending also on the mass of the satellite, the migration timescale is from 10 8 to 10 5 times slower than the orbital period, as shown in Section 4.1). As a result, since the orbital velocity of gas is sub-Keplerian (Pringle 1981), then any effect of azimuthal variations of surface density and aspect ratio on migration are averaged in time and negligible in the model. Unfortunately, this does not apply to dust, as its orbital velocity is considered to be Keplerian, but in Section 4.4 we show that the accretion process is also azimuthally averaged, this time by the differential velocity of the dust within the feeding zone of a satellite.
The dust-to-gas ratio ( ) of CPDs is still unknown. Therefore, we choose the dust-to-gas ratio of the disc as one of the initial parameters that are varied between the different runs of population synthesis (see in Section 3).
The initial disc density and temperature profiles are shown in Figure 1 and 2, respectively. With these profiles, the total mass of the CPD is found to be 9 × 10 −3 , and, in the outer part, the thickness of the disc is about one order of magnitude bigger than the one usually assumed for protoplanetary discs (∼ 0.5 versus 0.05).

Disc Evolution
The initial density and temperature profiles of the CPD are evolving in time in our semi-analytical runs, because the disc is dissipating and cooling. As in Cilibrasi et al. (2018), the gas density and temperature evolution is treated as an exponential law where the dispersion time is taken to be disp = 10 5 years. That means that the disc is completely dissipated by the time the simulations end, i.e. after 10 disp = 10 6 years. Consequently, the model is able to capture the system evolution even after the gas is completely dissipated (i.e. in the debris disc phase). The minimum background temperature at Jupiter location is assumed to be min = 130 (as in Miguel & Ida 2016) and the cooling time is calculated via radiative cooling (see Wilkins & Clarke 2012;Cilibrasi et al. 2018) with a timescale of cool 1.6 × 10 5 years.
The dust profile evolution is more complex. Overall, the dust profile is dissipating with the same time evolution as the gas component, but the interaction with satellites causes some other effects that are taken into account. First of all, as described in Section 2.3.2, the satellites accrete mass from the dust component, with a certain accretion rate within their feeding zone. Those cells that are part of the feeding zone will experience a density decrease that we call Σ acc . Furthermore, the same dust refilling mechanism used in Cilibrasi et al. (2018) to model the solid accretion rate from the PPD to the CPD has been implemented, i.e. the dust coming from the PPD will replenish the gaps in the CPD created by solid accretion of the seeds in a timescale ( ref ). This effect tends to bring the dust density to an equilibrium solution, that is decreasing exponentially. This timescale, as the dust-to-gas ratio of the disc, is an unknown parameter and will be treated as a free parameter in this model.
Mathematically, the equation for dust evolution that is solved in the model is the following: where the first term represents the exponential evolution, the second is the source term due to accretion, and the third term represents the refilling mechanism. With no accretion, the solution reduces to Σ( , ) = Σ 0 ( ) − / disp .

N-body interactions
In our model, protosatellites are treated as point masses orbiting around the central planet, interacting with the CPD and with each other in 3D, using a 3D N-body integrator (Grimm & Stadel 2014).
To begin with, each simulation starts with 10 to 20 embryos randomly distributed in the disc between 30 and 200 , and the initial semi-major axis distribution is such that log 10 is uniform. All embryos initially have 0 eccentricity and a random small inclination, distributed linearly between −0.01 and 0.01, i.e. the typical values for the current Galilean satellites. Any lost satellite (because of collisions or ejections, as explained further in this section) is replaced with a new embryo in another random position in the disc, making the initial number of embryos also the maximum number of particles in the disc at a given moment. The initial mass of satellites is chosen to be 10 −7 and a new embryo is created only if the total mass of the dust present in the disc in that moment is at least 10 times the mass of the embryo itself.
Within the N-body integrator, satellites are considered as particles orbiting a central fixed particle with = 1 . As the central particle is fixed, we consider it as an external potential, and all the sums in this Section will refer to the N satellites only. This simplifies the algorithm, since we do not need to take into account the effect of the planet movement in the disc, and it is reasonable, since the mass of the satellites is always significantly lower than the planetary mass (see in Section 3.1 for the total mass distribution).
The Hamiltonian of the system is then where is the momentum of the satellite , is its mass, its distance to the centre and is the distance between the particles and . This Hamiltonian is integrated following the algorithm shown by Saha & Tremaine (1994). First, a leap-frog time-integrator is implemented, dividing the Hamiltonian into the Keplerian part (the drift) and the Interaction part (the kick), i.e.
This is more accurate than dividing the Hamiltonian into kinetic and potential energy, because the Interaction part is almost always a small correction to the Keplerian motion and this gives a more accurate method with the same cost. Analytic solutions for the Keplerian motion are well known and this model uses the same and function method used in GENGA (Grimm & Stadel 2014) and described in Danby (1988).

Individual time-steps
In the model, the required precision is that any satellite needs at least 25 time-steps to complete one orbit. If the time-step were the same for all satellites, that would have been computationally too expensive due to the different orbital velocities of the inner and outer satellites within the disc. Instead, we split the Hamiltonian into multiple parts, each one related to one satellite This allows to solve the different parts of the Hamiltonian individually with different time-steps following the algorithm in Equation (14) in Saha & Tremaine (1994), thanks to the properties of the operators defined above: The time-steps have to be chosen so that they are all proportional to a power of 2, i.e. Δ = 2 Δ 0 , where Δ 0 is chosen to be the orbital time of the inner orbit ( = 1 ) divided by 25, i.e. the required precision.
is chosen in order to have Δ as the minimum time-step with the required precision, i.e. less than the orbital time of the ℎ satellite divided by 25. In practice, if is the period of the ℎ particle, then = int log 2 T i 25Δt 0 . During the evolution of the system, satellites may vary their time-step because they move through the disc. In order to keep the code stable, a single satellite can only decrease its own time-step when migrating. Finally, the algorithm provided by Saha & Tremaine (1994) needs the satellites to be ordered from the smallest to the larger Δ . This requires satellites to change also their index order in the code when migrating and changing their individual time-step.

Close-encounters
The algorithm works well as long as is a small correction of . If the two become comparable, then the error of the algorithm, as given by the Baker-Campbell-Hausdorff formula (Saha & Tremaine 1992, 1994, becomes one order of magnitude larger, and that happens when close encounters between satellites occur. This is managed implementing a close-encounter treatment from Chambers (1999), following the example of GENGA (Grimm & Stadel 2014). In case of close-encounters the Hamiltonian portions change to where is a changeover function calculated for each pair of satellites. It has to be 1 when satellites are far away and it needs to go to 0 when they are coming closer. In particular where crit = 3 Hill, and Hill, = as suggested by Duncan et al. (1998). If close-encounters were involving also small satellites, a correction depending on velocities needs to be added for calculating crit , as explained in Chambers (1999), but that is not the case in this model, as described below. The shape of is chosen to be continuous and derivable twice in the whole domain. Furthermore, we automatically assign = 1 when at least one of the two satellites or has < threshold , where the threshold mass is chosen to be 5 × 10 −6 . This is because close-encounters and collisions between small embryos are not significantly influencing the outcomes of this model. When two or more satellites have < 1, the algorithm needs to change because { , } ≠ 0 and the Keplerian parts of the Hamiltonian cannot be solved individually and analytically. The Hamiltonian where , are the indexes of satellites involved in the close-encounter, (i.e. the ones for which < 1) can be solved, for example, implementing another leap-frog algorithm with a smaller time-step. As in Grimm & Stadel (2014), we decided to implement the Hermite Integrator with Ahmad-Cohen Scheme presented in Makino & Aarseth (1992). This method is fourth-order because it computes the force and its derivative analytically, and then construct a third order polynomial interpolation. The method needs two parameters to compute the appropriate time-step at the beginning and during the integration. We chose these parameters to be = 0.02 and = 0.01, as defined in Makino & Aarseth (1992).
This close encounter treatment, as presented in Chambers (1999) and as used in Grimm & Stadel (2014), needs to be implemented in a democratic coordinate system, which consists of heliocentric (in our case planetocentric) positions and barycentric velocities. Since we have a central fixed particle (with infinite inertial mass) the democratic coordinates coincide with the Cartesian coordinate. Combining together individual time-steps and close-encounters is fundamental for this model in terms of performance and accuracy, but a synchronization problem may occur if two satellites that come closer do not have the same time-step, because it is not possible to integrate close-enc . If the two particles are synchronized, i.e. their internal clock is the same (see Saha & Tremaine 1994 for details), there is no such issue. Then, the longer time-step is changed to the shorter timestep and the calculation can go on as described above. If that is no longer the case, then the code comes back to the last moment where the two satellites were synchronized, changes both the time-steps and starts over from there. In order to do that, the code is keeping in memory some back-up states of the system long enough to deal with these situations.

Collisions
During the N-body integration, collisions are also possible. Assuming a density around the mean density of Galilean satellites ( 2.5 / 3 ), a radius is assigned to each satellite. When the distance between two satellites is less than the sum of their radii, then they are considered colliding. Then, the two satellites are replaced by one, having the sum of the individual colliding moons' masses, situated in their common centre of mass and proceeding with the sum of the two momenta. In other words, the collisions are perfectly inelastic, meaning that the momentum is conserved, but the energy is not. In this model, we did not consider the possibility that satellites could also experience disruptive collisions, fragmentation or only partial merging. That would require a significantly more detailed modelling resolving their radii, which is beyond the scope of this paper.
Satellites can also collide with the central planet (i.e. get engulfed by it). That happens when the radial distance of a satellite from the centre is less than 1 . In this case the satellite is deleted from the simulation and its mass is considered as lost into the planet. The third possibility is that a satellite assumes an eccentricity that is ≥ 1. In this case the orbit is parabolic or hyperbolic, hence ejected from the system.
In the above three scenarios, when one satellite is lost, it is replaced by another embryo, randomly located in the disc following the same initial prescriptions regarding the initial position and mass, as described before.

Migration and Damping
While interacting with each other within the N-body integrator, satellites also interact with the disc itself. This interaction causes migration, i.e. a change in the semi-major axis of the orbit (usually inward), and eccentricity and inclination damping. In this model, these effects are implemented following the prescriptions by Cresswell & Nelson (2008). These additional accelerations are included in the Interaction (Kick) part of the Hamiltonian as where mig , ecc , and inc are functions of eccentricity and inclination defined in Cresswell & Nelson (2008), and wave = Type II migration replaces type I when < 1 and the migration parameter is calculated following the Pringle (1981) model and the damping prescribed by Syer & Clarke (1995) More in detail, the two regimes are connected using the junction function proposed by Dittkrist et al. (2014) where ( ) = (1 + 30 ) −1 .

Accretion
The satellites accrete dust from the CPD, acquiring mass this way too. In particular, satellites that are moving across the disc will always accrete the solids from within their feeding zone. Following Greenberg et al. (1991), the solid accretion rate is: and the feeding zone extends up to 2.3 Hill . When a satellite accretes, the cells within its feeding zone are going to have lower dust density as a result. The same way, accretion happens as long as the cells within the feeding zone have enough mass to be accreted. The cells are then replenished by the dust refilling mechanism explained in Section 2.1.2.

RESULTS
The satellite formation model described in Section 2 was run following a population synthesis approach (Ida & Lin 2004;Benz et al. 2014;Mordasini 2018). Three parameters were varied randomly with different values for each simulation, producing different initial conditions and therefore, different outcome. The three parameters were the dust-to gas ratio, the dust refilling timescale ref and the number of initial satellite seeds. The idea of populations synthesis is to statistically analyse the outcome of random initial parameters on the resulting population. In this case the ranges for the parameters have been chosen to be The distribution shapes for these parameters were chosen in different ways. The initial number of seeds has been chosen uniformly in the range between 10 and 20. The other two parameters were distributed finding a compromise between the same logarithmic distributions used in Cilibrasi et al. (2018) and the computational time. In fact, systems with a higher dust-to-gas ratio and a faster refilling, have more massive satellites and, consequently, the close-encounter integrator is much slower. The resulting distributions are shown in Figure  3. A total number of 2309 systems were simulated, enough to reach convergence on the results. We tested higher simulation numbers too, but the shape and characteristic parameters (mean, variance) of the distributions did not change anymore above 2000 runs. Once all the simulations have been completed, the idea of a population synthesis is to look at the statistics of the results, gaining a general understanding of how satellite systems around Jupiter-like planets might look like. On the other hand, since this analysis depends on the initial parameter distribution, another approach is looking at the dependence of the results on these parameters, allowing pattern identification and a deeper understanding of the significance of the results.

Retained and lost mass
During a single simulation, many satellites were engulfed by the central planet, collide with each other or get ejected from the system. First, we looked at the final mass of single survived satellites. From the histograms and the cumulative distribution in the first plot in Figure 4, it is clear that small satellites are dominant and that the mass distribution is strongly peaked towards lower values. Nevertheless, most of this light satellites are just seeds that are produced throughout the simulation and have little effect on the real architecture of a satellite system. For example, the second plot of Figure 4 shows again the masses of single survived satellites, but the histograms are built weighting the satellites with their mass. In practice, the value of each column represents how much mass over the total mass of all the survived satellites is contained in that particular bin. In this case, we see that massive satellites are now dominant, meaning that, even though numerous, small satellitesimals are not important in determining and characterizing the structure of the survived systems.
Consequently, we identified the total mass of survived satellites as a more significant measure. Figure 5 shows that the low mass systems are more likely than higher mass systems. If we consider the Galilean system mass as a threshold (red line in Figure 5), i.e. 2 × 10 −4 = 0.06 ⊕ , then 85% of the produced systems are less massive than this threshold, while 15% of the systems are more massive. Figure 6 shows the total satellite mass distribution within the resulting moon-systems, but ignoring satellites less massive than the chosen thresholds. The figure shows that if the contribution of all the small satellitesimals, the shape of the distribution will change. In particular, the most massive systems do not vary a lot, because small satellites give only a small contribution, while the distribution changes quite much in the left tail. This is because low-mass satellites have a bigger impact on less massive systems. In particular, if only Galilean-like satellites ( ≥ 10 −5 ) are considered, the peak of the (green) distribution is very close to the total Galilean mass (2 × 10 −4 ). On Figure 7 the total mass of satellites that were engulfed by the planet is shown. These moons contribute to the heavy element content of the gas giant planet envelope, and some probably sediment down to the core, possibly adding fuzziness to the core boundaries. This "metal-pollution" happens only in 10% of the cases. In 10% of them, i.e. a total of 1% of the total cases, the polluting solid-mass is In the second case, the distribution has been weighted by the mass of satellites itself, so that the sum of the values of all the bins equals 1. This highlights that massive satellites are lower in number, but more important in terms of mass budget. comparable to 1 ⊕ . From this we can conclude that solid accretion onto the planet is not as efficient as in our previous work Cilibrasi et al. (2018). The difference rises due to the absence of a dust trap in this work (Section 4.2) and the 3D treatment of close-encounters and collisions (Section 2.2) instead of 1D, that leads to less massive moons and, consequently, slower migration.

Number and locations
The number of survived satellites and their semi-major axes were also examined. Figure 8 shows the probability distributions of the number of survived satellites above certain threshold masses (10 −7 , 10 −6 , 10 −5 ). These thresholds are used to distinguish again satellite seeds that are always produced by the model, from fullgrown moons, given the lack of a proper definition between the two mass categories. The distribution in Figure 8 with ℎ = 10 −7 (i.e. all satellites and satellite seeds) are for comparison purposes only. As the mass threshold is increased from the base value ( ℎ = 10 −7 ), there is a greater number of systems that cannot build any satellites above the threshold (20% for th = 10 −6 ,  60% for th = 10 −5 ), and the shape of the distribution changes as well. For example, the case with the threshold of th = 10 −5 (i.e. counting only moons that are on the individual Galilean moon mass regime), about 18% of the cases have 3, 4, or 5 satellites, in agreement with the Galilean system. This is trivial, since the same CPD can only produce smaller amount of heavier moons, or more numerous low-mass moons.
In order to reproduce a Galilean-like system, it is also necessary to check if the satellites migrate from the formation location to the typical Galilean locations. The distribution of semi-major axes of satellites from our population (in Figure 9) shows that the distribution is mostly concentrated around the formation location because light satellites do not migrate much. When increasing the mass-threshold that we examine, the distribution tends to spread in the semi-major axis range a bit, especially towards the planet location. Satellites in the current Galilean-moons' locations increase from 15% to 20%.  Furthermore, in Figure 9 there is a peak in the distribution in the outer part of the seed inserting region. This happens because migration is less effective in the outer part of the disc (the timescale is longer if computed as in Section 2.3) and that causes satellites to migrate less when they form in that outer region. This is also the reason why the distribution spreads more towards the inner disc than towards the outer part when changing the mass threshold. The semi-major axes of satellites are also linked to the resonances that may form in a satellite system, similarly to planetary systems. In fact, a mean motion resonance (MMR) occurs when orbiting particles or bodies exert a periodic and regular gravitational influence on each other, and that happens usually when their orbital periods, that is given by their semi-major axes, are related by a ratio of small integers. In this analysis, an approximate treatment has been implemented, and a pair of consecutive satellites and is considered as in a resonant configuration : , where , ∈ , or in this case , ∈ {1, 2, 3, 4, 5}, when their period ratio is closer than 0.5% to the ratio  / (Hands et al. 2014). We perform the analysis for all the satellites (about 24 thousand pairs) and for all satellites with > 5 × ℎ = 10 −6 (about 5 thousand pairs), in order not to take into account the effect of small satellitesimals. As shown in Figure 10, when all fullgrown satellites and satellite-seeds are included in the examination, 6.3% of the pairs were found to be resonant, while considering only the larger moons, there are 8.5% of resonant pairs. In particular, 1 : 2 and 2 : 3 are the most common configurations.

Architecture of systems
So far, distributions for moon-masses, -numbers and semi-major axes have been examined whether their values are close to the Galilean system. The previous analysis has two limitations: it depends on the distributions of initial and boundary parameters (dust-to-gas ratio, dust refilling timescale, initial number of satellite seeds), secondly, it does not quantitatively tell how similar a single produced system to the Galilean one. In order to overcome these limitations, it was defined how "close" a system is to the others or to the Galilean system, following the approach of Alibert (2019). In particular, Alibert (2019) shows how a mathematical distance between the similarity of systems can be defined. First of all, each moon-system is assigned a 2D function depending on two variables, a mass and a semi-major axis , defined as with Here is a function calculated for each satellite ( ) in the ℎ system, and are the mass and the semi-major axis of the satellite , and are two arbitrary parameters that have to be tuned in the analysis, and and are the two aforementioned functional variables. The functions are defined over a reasonable domain, i.e. (10 −7 ≤ ≤ 10 −2 , 1 ≤ ≤ 1000 ), in order to cover all the possible combinations.
In practice, the 's are sums of 2D Gaussian profiles, one for each satellite in the system, with and defining the width of the Gaussians in the two directions. ( ) is an arbitrary weight function that controls the height of the Gaussian peaks. In fact, the goal is to compare systems giving more importance to massive satellites, and that is why more massive satellites should have higher Gaussian profiles. We identified two possible weighting functions, depending on the desired application: log ( ) = log 10 ( / ) + 7 with the first function being the one used in this analysis. Once every system has been associated to a function ( , ), a norm of these functions has been defined as || || = ∫ ∫ | ( , )| 2 d log 10 M d log 10 a and as a consequence we have a distance, in arbitrary units, This is proven to be an actual distance, fulfilling all the necessary conditions, as the functions are part of an 2 space. Furthermore, we found convenient in the analysis to normalize the 's in order to have || || = 1 ∀ . The first step of the comparison procedure is to check the distance of the individual moon-systems that came out from our population synthesis to the Galilean one. We defined a for the Jovian case and computed the distance from all other systems , using = 0.1 and = 0.1. This means that two moons are considered very similar (i.e. their Gaussian profiles almost coincide and the difference gives a quasi-null integral) when they differ up to 1/10 of order of magnitude in mass and semi-major axis. All the distances are shown as a function of dust-to-gas ratio and dust refilling timescale in Figure 11, as the initial number of seeds turns out not to have an effect, because the number of massive satellites is determined by the amount of solids in the CPD. Given that normalized 's can only have a distance between 0 and √ 2, there are high distance points in the upper left corner and in the lower right in Figure 11, while there is a band, that is highlighted in red, where the distances from the Jovian system get smaller. In particular, measuring the slope of the band, we identify Figure 11. Distance of all the systems from the Galilean satellites as a function of dust-to-gas ratio and dust refilling timescale. The unit of distances is arbitrary in this analysis framework. The red band highlights Galilean-like systems.
with a peak around ref 2.2 = 2 × 10 8 yrs. We can conclude then that the fundamental parameter, used to predict whether a system will produce Galilean-like satellite or not, is actually the ratio Γ ≡ ref 2.2 , as also shown in Figure 12. This relation can be verified visually from the plots, but it can also be derived more formally with the socalled Randomized Dependence Coefficient (Lopez-Paz et al. 2013), as explained in Appendix A.
In order to verify the aforementioned relation again, and in order to understand how really similar to Galilean satellites the produced systems are, we performed further analysis. In fact, being the distance in arbitrary unit by definition in this analysis framework, the previous analysis, shown in Figures 11 and 12, tells how to pro- Figure 13. Results of the t-SNE algorithm (see details in the text). The systems are reproduced in a 2D plot. The axes do not have any physical meaning, as the t-SNE approach is just a way to visualize high dimensional data in a 2D space, preserving local structures. Similar systems appear close and different systems appear far away. The green cross represents the Galilean system. duce the most Galilean-like systems in the model, in terms of their properties, but it is not clear how similar they are. It is necessary to understand whether, in those conditions, the distance from the Galilean system is similar (or even less) than the typical distance between the moon-systems produced by our population synthesis. Here, we are working in the infinite dimension space of 2 functions in the chosen domain, hence it is not possible to easily visualize distances. Following the Machine Learning approach of Alibert (2019), we also apply a T-distributed Stochastic Neighbor Embedding, or T-SNE, algorithm (Van der Maaten & Hinton 2008; Van der Maaten 2014) to our data. The T-SNE is a nonlinear dimensionality reduction technique meant for treating high-dimensional data for visualization in a low-dimensional space, typically a 2D arbitrary-unit space. In particular, minimizing a cost function, the algorithm links each highdimensional object with a point in the 2D space so that similar objects are linked to close points and dissimilar objects are linked to distant points with higher probability. This procedure prevents from the typical mistakes of projections, where distant objects could be projected as being close to each other in a 2D plane.
The t-SNE algorithm was applied to our data using the python package Scikit-Learn (Pedregosa et al. 2011), choosing a perplexity = 50, which is a parameter related to the number of nearest neighbors and balances the sensitivity of the algorithm to the local and global features of the data-set. The result is shown in Figure 13, in which the systems are reproduced in a 2D plot (the two axes do not have any physical meaning), so that similar systems appear close and different systems appear far away. The colour represents the value of Γ = ref 2.2 . First, the plot shows that the Galilean system (green cross) is well embedded in the cluster of systems. This means that the Jovian system is actually reproducible by our model in certain conditions. That would not be the case if the cross was out of the pipe-shaped cluster. Figure 13 shows also that the parameter Γ is actually controlling the architecture of the final outcomes of the model. In fact, following the cluster from one end to the other, we have a continuous and gradual change of that parameter, with the Galilean system being close to Γ = ref 2.2 = 2 × 10 8 , as predicted before.

Irregular satellites
Irregular satellites in the Solar System have large inclination and eccentricity and thought to be captured by the planets (instead of being formed around it). Our population synthesis also produced high inclination and high eccentricity moons, hence we examined how frequently they are formed within our framework. In our formation model, the moons are forming in a CPD around the planet, but the dynamical interactions can toss them to highly inclined, elliptical orbits. We therefore checked whether some of the irregular satellites of the Solar System can be reproduced from traditional moon-formation processes. In our analysis, we consider a satellite as irregular when the eccentricity is ≥ 0.1 and the inclination is sin( ) ≥ 0.1, following the definition presented by Jewitt & Haghighipour (2007). Irregular satellites around Jupiter and Saturn also have also values for eccentricities ≥ 0.1 and inclinations ≥ 0.1. With this definition, in the outcome of our population synthesis we found about 1% irregular satellites (Figure 14). In particular, our 2309 runs produced 26293 satellites, of which 262 are irregular. Among these, only one retrograde satellite was found, the others were prograde. On the other hand, among the Solar System irregular satellites, the retrogrades are far more common than the progrades. Figure 14 shows also the masses of the irregular satellites with a colour scale. The produced irregular satellites in our framework can be quite massive. In particular, about 4% of the total satellite mass is belonging to the irregular satellites. The mass histogram in Figure 15 shows that the peak of the distribution is at about 10 −5 , that is on the order of the Galilean satellites' mass. This is quite different if compared to irregular satellites in the Solar System. In fact, Jupiter, Saturn and Uranus (Neptune is believed to have had a more complex history in the satellite formation process, as shown by Agnor & Hamilton 2006) usually have a lot of irregular satellites in number (90%, 70% and 33%, respectively), but carrying a very small fraction of the total mass (0.001%, 0.007%, 0.04%, respectively).
In conclusion, our traditional moon-formation scenario within a CPD can reproduce only a small fraction (≤ 1%) of the irregular satellites known in the Solar System, highlighting yet again the capture origin of irregulars.

Observability of exomoon-systems
Since our framework produces a population of satellites that formed in CPDs similar to Jupiter's disc, it is possible to check what fraction of the produced moons would be detectable with current and nearfuture instruments. Here we consider only the required sensitivies, not the observational time related feasibility. The planet in our model is at 5.2 AU from its Solar-equivalent star (Jupiter's distance), for which having 3 transits measured would require ∼ 36 years.
For each system we computed which sensitivity is needed in order to have more than 50% probability of observing that moon-system (assuming that the planet is transiting in front of the star from our point-of-view). The reason for this 50% probability threshold is because the moons are not detectable during their entire orbits around the planet, due to the trivial geometrical effects. The sensitivity is considered as the instrumental capability to distinguish a transit curve from the noise. For example, an instrument is considered to have a sensitivity of 10 −5 if it is able to spot a 10 −5 single-transit dip in a stellar light curve. In this analysis, we identify how many moonsystems it is possible to detect given a sensitivity-limit, i.e. it is possible to infer a detection probability around a Jupiter-like planet as a function of the sensitivity. The details about the computation are in Appendix B. Figure 16 shows the distribution of the minimal sensitivity that each system needs in order to be detected. The black curve on the figure represents the detection probability over the whole population with the given sensitivity. As expected, this is 1 for very low sensitivities, because it allows to detect all the satellite-systems in our model, and it goes to 0 for higher sensitivities. A sensitivity of 10 −5 would allow to detect about 14% of the systems produced within our model, while we find the median value at 2 × 10 −6 (a sensitivity of 10 −6 would allow to detect more than 60% of the moon systems).
In order to have a more general insight into the detectability of moons, we also investigate the magnitude of the TTV (Transit Timing Variation) that the satellites produced by our model would have on a Jupiter-like planet transits. This effect has been calculated using the formula by Kipping (2009) and assuming zero eccentricity, i.e.
where is the semi-major axis of the satellite, the one of the Figure 16. Distribution of sensitivities needed to detect a moon-system. The black curve represents the fraction of systems that we would be able to detect with the given sensitivity. The analysis is on all 2309 systems. planet, and the masses of the satellite and the planet, the orbital period of the planet.
The TTVs calculated for our population are shown in Figure 17. The distribution spreads from 0.01 seconds to over 100 seconds, has a mean of 2.56 sec and a median of 0.26 sec. Furthermore, as the cumulative distribution in Figure 17 shows, about 6.2% of the satellites would produce a TTV greater than 10 sec.

DISCUSSION AND CAVEATS
This Section describes the major caveats and biases of our framework, identifying their effects on the results and comparing them to previous literature. In particular, we discuss the choices of the migration formulas and the chosen dust density profile, we compare our model with a simplified 1D model, we analyse the differences between a Figure 18. Migration rates as calculated with the formula from Jiménez & Masset (2017). The rate is the inverse of the migration timescale multiplied by the orbital frequency, i.e. it is the orbital timescale divided by the migration timescale. Figure 19. Migration rate comparison between formulae from Paardekooper et al. (2010Paardekooper et al. ( , 2011 and Jiménez & Masset (2017). The rate is again the inverse of the migration timescale multiplied by the orbital frequency, i.e. it is the orbital timescale divided by the migration timescale.
1D and a 2D disc model, and we identify future developments of the model.

Migration prescriptions
As explained in Section 2.3, the parameter for type I migration has been taken from Jiménez & Masset (2017), because it was derived from 3D non-isothermal simulations, and can be applied to lowand intermediate-mass planets in optically thick discs. On the other hand, other migration formulae can be found in the literature (Tanaka et al. 2002, D'Angelo & Lubow 2010, Paardekooper et al. 2010, 2011, with the latter being the one used in our previous work, in Cilibrasi et al. (2018). The Paardekooper-formula was derived from simulations with low-mass planets embedded in 2D non-isothermal discs, and here we compare it with the one by Jiménez & Masset (2017), showing that they give very similar migration rates in our domain. Therefore the differences between this work and Cilibrasi et al. (2018) are not due to the different migration timescales. The migration rate, as calculated following Jiménez & Masset (2017), is shown in Figure 18 as a function of the satellite semimajor axis and mass. The plot shows that migration is always inward (always negative) and it is always very slow. In fact, the ratio between the orbital timescale and the migration timescale is between 10 −9 and 10 −5 . Figure 19 shows a comparison between the migration rates as calculated with Jiménez & Masset (2017) and Paardekooper et al. (2010Paardekooper et al. ( , 2011. In this case we see that the ratio of the two rates is always around 1, with a maximum spreading between 0.25 and 4, in very localized areas in the parameter space. This means that the two rates are always comparable well within one order of magnitude and we do not lead to significant differences for our population synthesis outcome. The plots also show the transition between the two different slopes of the temperature profile.

Formation time and dust-trap effect
In the resulting moon population of our model, the formation timescales of satellites were also examined. Because the different massive moons grown on different timescales, we chose to look into how much time it takes to grow to Europa's mass, which is the lightest Galilean moon (2.53 × 10 −5 = 8 × 10 −3 ⊕ ). The formation timescale histogram in Figure 20 shows 1792 satellites, that have reached this threshold mass. The histogram shows a one-peak distribution with a long tail on the left, toward smaller timescales, and a shorter tail on the right, given by the few cases in which satellites grow due to collision until the end of the simulation (after 10 6 years). Furthermore, the histogram can be fitted with a Cauchy distribution ∝ 1 + log 10 t form − 2 −1 with = 5.08 ± 0.01 and = 0.37 ± 0.02, being the peak of the distribution at about 10 5 years. In particular, about 53% of the satellites take more than 10 5 years to reach the Europa's mass.
Hundred-thousand years is a special threshold for satellite formation timescales, given the internal structure of the Galilean moons. In particular, we know that the first three of them (Io, Europa, Ganymede) are differentiated, while Callisto is not. These different features are thought to be linked with the formation timescale. A fast formation and accretion can cause a satellite to melt, because of the high rate of energy absorption, allowing heavier elements to sink and form layers. On the other hand, a slow formation process would not cause melting, preventing the moon from differentiation. The threshold accretion time that divides the two scenarios is thought to be 10 5 years (Stevenson et al. 1986;Canup & Ward 2002).
Compared to our previous work, Cilibrasi et al. (2018), the peak of the distribution in Figure 20 is considerably shifted to the right, showing much slower formation timescales this time. This is mainly due to one difference in the two models, i.e. the dust distribution in the CPD (other minor differences, such as the 3D vs 1D dynamical integration, turns out to have a minor effect and are investigated in Section 4.3). While the model presented in this paper consider the dust distribution to have the same shape and slope of the gas distribution, the model in Cilibrasi et al. (2018) used a dust distribution derived from a 1D dust evolution model by , that shows a peak in the distribution at about 85 from the planet. This location is a dust trap and was considered to be the main location where seeds formed within the CPD. Here, the moonlets acquired most of their mass, which happened very quickly, given the high concentration of dust. In this paper, the smoother dust density profile allows seeds to form everywhere in the disc, but that leads to slower growth rate. Consequently, this difference also explains why we there is no frequent sequential formation in this work. Slowly growing satellites do not experience fast Type I migration, and hence less of them gets engulfed by the planet.

1D vs 3D
In order to understand the importance of using an N-body integrator in our model, we compared our results with runs performed in a fully 1D population synthesis model. In this test, the disc evolution was treated exactly the same way, while the satellite dynamics was treated as in Cilibrasi et al. (2018), assuming circular orbits, migration rates by Jiménez & Masset (2017) and resonant trapping as described by Ida & Lin (2010). Collisions were triggered when two satellites approached one another at a distance smaller than the sum of their respective Hill-radii. Because the computational time is orders of magnitude lower than of the 3D runs, we were able to perform a 1D simulation starting from the same initial conditions and parameters for each of the 3D simulations.
We compared the distributions of the total satellite mass in each system ( Figure 21). The difference is that the satellites can grow to larger masses in the 3D model relative to the 1D model. The reason for this is mainly how collisions are treated. In 1D the collisions are trivial, while in 3D the collisional cross-sections are significantly smaller. In the first case, this accelerates the mass growth of moonlets, which then increases their Hill-radius, leading to a further increase in the collision rate. However, because of the larger masses, satellites also migrate faster towards the planet, eventually getting engulfed by the gas-giant. On the other hand, in the 3D model, this runaway growth is not occurring, which in turn leads to a lesser role of migration in the orbital evolution. At the same time, in the 3D model, collisions can scatter small satelletesimals towards the planet, that would not be affected by inward migration given their low masses. Scattering is completely absent in the 1D model due to the dimensional limitation. Figure 22 shows the mass distribution of satellites that are lost into the central planet. This figure again demonstrates the different outcomes of the 3D and 1D models, similarly to the explanation above. The curve of the 3D model shows two peaks, one for small satellites, which are delivered to the centre by scattering, and one for larger masses, which are the result of inward migration. On the Figure 21. Distribution of the total mass of survived satellites (at the end of the simulations) in each system in the 1D and 3D model. The distribution is normalized so that the integral is 1, as we consider the same number of systems in the two cases. other hand, the curve for the 1D model shows only one peak, which is caused by migration only, since there is no inward scattering in 1D.
In order to understand more qualitatively the difference between the 1D and 3D treatments, the same machine learning framework was used as described in Section 3.3. Thanks to the fast running time of the 1D simulations, we performed 20 simulations for 300 different combinations of , ref and init , each time starting from an initial random distribution of seeds. This way it is possible to estimate the spreading of systems simulated with the same parameters, i.e. the dispersion of results that is naturally caused by the randomness of the code (initial seed positions), and compare it to the distance between the 1D and the 3D model. Using Γ = ref / 2.2 as the reference parameter, we show this comparison in Figure 23. The Figure 23. The blue dots represent the distance between the 1D and the 3D results in each system, while the red area is the distance dispersion range of systems generated by the same initial parameters in the 1D model distance between 1D and 3D results (blue dots) is compared with the typical dispersion range of distances between systems generated from the same parameters (red range), making clear that the two values are comparable only in the range where ref / 2.2 is high, i.e. where satellites grow and migrate less, while they are well-distinguished for lower values. This means that the 1D and the 3D models produce significant differences.

1D vs 2D Accretion
We also tested what differences it would make if the CPD was treated in 2D, instead of 1D. Since our disc model is coming from a 3D thermo-hydrodynamical simulation, it is possible to consider the midplane of this disc, and evolve it keeping its 2D structure (e.g. the spiral arms present in the disc).
We verified in the 2D case too, that migration timescales are much longer than orbital timescales. This means that assuming a 1D profile is giving just a negligible correction to migration and damping calculations compared to a 2D model, since the orbital motion of the Keplerian satellites through the sub-Keplerian gas is already averaging all disc quantities azimuthally.
On the other hand, solid (dust) accretion is different in 2D and 1D. In fact the satellites and the solids are in principle orbiting with the same Keplerian speed and do not necessarily show the same azimuthal averaging process. This means that a moonlet does not have access to all the dust in its orbit at once, as in the 1D model, but only to the amount in its 2D feeding zone ( feed = 2.3 Hill ; Greenberg et al. 1991). Nevertheless, one can expect the differential velocity of all the dust in the feeding zone to produce some azimuthal average as well. In order to understand whether or not this average process actually happens, we need to calculate how different accretion is in the two scenarios (1D versus 2D). The process is not linear, since a satellite would accrete more if it is more massive (since the Hill radius would be larger, hence its feeding zone as well). To test the difference that this makes, we set up 10'000 pairs of simulations. Each pair had the same satellite with a mass randomly chosen between 10 −7 and 10 −3 , a semi-major axis between 10 and 300 , embedded in a disc with a dust-to-gas ratio between 0.001 and 0.5 and refilling timescale between 10 2 yrs and 10 6 yrs. In each of the two simulations in the pair, the satellite was let free to accrete mass from the dust of the disc, until one of the two following conditions happened: • the time in the simulation had reached 100'000 orbits at the location of the satellite • the satellite had migrated more than 0.1% of its semi-major axis because of migration, calculated as Section 2.3.
In the pairs, one simulation was performed with a 2D orbiting (Keplerian) dust grid and the other simulation was performed with a 1D grid, with the results being compared afterwards. In the analysis, the initial mass of the embryos and the dust-to-gas ratio of the disc turned out to have the biggest influence, while the semi-major axis and the dust refilling timescale did not have a fundamental effect.
The dependence on the initial mass of the satellite seeds and the dust-to-gas ratio is shown in Figure 24. As expected, higher dustto-gas ratios produce higher differences in the final mass ratio: the more massive moons accrete faster and create a bigger discrepancy. On the other hand, embryos have a different evolution in the very first phases, when the mass is still small. Once satellites become more massive, then this effect is much smaller. The plot shows that the difference between the 1D moon masses (being the larger) and the 2D ones is never more than 15%. In the mass ratio histogram in Figure 25. Distribution of the relative difference between 1D-and 2Dcalculated moon-masses. The peak is at about 0.2−0.5% while the maximum difference is at about 15%. The two colours identify satellites that stop the accretion because of migration and the ones that stop because of reaching N=100'000 orbits (i.e. end of the simulation). Figure 25 between the 1D and 2D cases, we see that the difference is more than 1% in 19% of the cases, while it is between 10 − 15% in 2% of the cases. This means that having a 1D dust profile instead of a 2D would produce a significant difference in dust accretion only in 2% of the cases. Once the small satellites reach higher masses and start to migrate through the disc, then the difference becomes even smaller.
On Figure 25 the high peak in the mass ratio histogram (that is linked to the over-density areas in Figure 24) is due to a bias in the set-up of the simulations. Even with high dust-to-gas ratios, in cases when dust refilling is slow, the disc is able to provide a maximum of about 10 −4 of dust-mass to the satellite before migration starts. Satellites often reach that maximum mass and either they stop their accretion (small initial mass), keeping the same Δ / until the end of the simulation, or the simulation stops because of the migration timescale constraint (big initial mass). This "saturation" value of Δ / turns out to be 5 × 10 −3 and it is mostly due to the exactly co-rotating dust grains.

Biases and Future Developments
In Section 2.3, the accretion model is a relatively simple one, but accretion itself indeed is a very important process, that affects the outcome. This can be seen from Figure 26, which represents the fraction of mass accreted from the dust disc over the total mass gained. It is clear, that most of the proto-satellites grow because of dust accretion, less because of collisions with other moonlets. The amount of dust accreted is of course affected by the density profile evolution, the disc dispersion time and the dust refilling mechanism (i.e. how much dust is getting into the circumplanetary disc from the meridional circulation). Instead of our oligarchic growth method, other accretion schemes could have been implemented too. For example, using pebble accretion (Ormel & Klahr 2010;Lambrechts & Johansen 2012), which was already applied in the satellite-formation models of Shibaike et al. (2019); Ronnet & Johansen (2020).
One of the main missing element, is the lack of dust-gas hydrodynamical simulations of circumplanetary discs, that can inform us Figure 26. The histogram shows the fraction of the total mass ( tot ) gained by dust accretion ( dust ). Most of the mass of moons gained by accretion from the dust disc, as opposed to collisions. The different colors represents the different moon-mass classes in the model. Small satellites (blue bars) grow mainly because of dust accretion (as shown by the peak around 0.9), while bigger satellites experience collisions more often, even though most of their mass still comes from dust accretion.
about the global and local dust-to-gas ratios and the dust flows based on particle size. Furthermore, our knowledge on the dust coagulation in the circumplanetary disc is also very limited . These effects of course very much affect any moon-formation model, since the satellites are forming from solids.
Because the base hydrodynamical simulations were boundary layer accretion simulations, there is no cavity between the planet and the circumplanetary disc in our model. Such a cavity would affect the formation of resonant chains (Ogihara & Ida 2009Sasaki et al. 2010;Fujii et al. 2017).

CONCLUSIONS
In this work, we investigated the formation and the evolution of satellites in a circumplanetary disc around a Jupiter-like planet, implementing a semi-analytical model for the disc evolution, and an N-body integrator for satellite-dynamics. We used a population synthesis approach, then statistically analysed the outcome.
In particular, the circumplanetary disc density and temperature profiles were taken from radiative 3D hydrodynamical simulations (Szulágyi 2017) and adapted into 1D and 2D grids. The dust profile was taken as a fixed fraction (dust-to-gas ratio) to the gas density, and the density and temperature profiles in the model evolve following an analytical (exponential) evolution. The proto-satellites formed and evolved interacting with the disc and consequently experiencing migration, eccentricity and inclination damping (Cresswell & Nelson 2008), and interacting with each other (using an N-body integrator with individual time-steps (Saha & Tremaine 1994) and a closeencounter treatment (Chambers 1999)). Satellites also accreted mass from the dust component of the disc, collide with each other, and migrate in the disc, while some of them get engulfed by planet. More than 2'000 individual simulations were run in a population synthesis framework, each one with a different dust-to-gas ratio, dust refilling timescale and number of initial seeds, chosen within appropriate limits (described in Sect. 2).
Our results show that the average satellite mass is below Galileanmasses, only 15% of moon-systems being more massive than the Galilean system. The moonlets grow slowly, and most frequently, it takes 10 5 years to reach Europa's mass. The migration rate is relatively slow, which means that only a few formed moons migrate all the way into the planet. In 10% of the cases moons are engulfed by the planet, but only 1% of the satellite-systems lose 1 Earth-mass or more into the planet. These lost moons and moonlets will contribute to the heavy element content of the giant planet's envelope. In all cases, the mass of the moons is mainly gained through accretion of dust, as oppose to collisions between moonlets. This highlights that the considered dust disc (in terms of dust-to-gas ratio, dust refilling timescale, etc.) might have a major role on the outcome of similar planet-formation models.
Checking for resonances, when all full-grown satellites and satellite-seeds are included in the analysis, 6.3% of them are found to be in resonant-pairs. When only massive satellites are considered, 8.5% of them are in resonance, with the most common configuration being (1:2) and (2:3).
If we consider only the Galilean-mass moons, 18% of our systems have 3, 4, or 5 satellites. Smaller satellites are of course more numerous. The eccentricies and inclinations of the resulting population was also examined. We defined an "irregular satellite" population (i.e. high eccentricities and inclinations) with having both of these parameters above 0.1. Only 1% of the total population belong to this irregular satellite-like group. Compared to irregular satellites in the Solar System, we found that ours are more massive and almost all cases (except one) prograde. This highlights that the moons formed in a circumplanetary disc cannot reproduce the irregulars at all, even if scattering is accounted for.
We used a machine learning approach to compare our population with the Galilean system, applying the t-SNE algorithm, similarly to Alibert (2019). This way, a "distance" between systems, meaning their similarity in the mass × semi-major axis space, has been defined and used to visualize and classify the outcomes. This approach demonstrates that the Galilean system is one possible outcome of the model and that the outcomes of the model depend almost solely on the parameter (Γ = ref / 2.2 ). In particular, a Galilean-like system is produced when Γ = 2 × 10 8 yr.
Our population of moons can be used to infer a potential exomoon population as well, around Jupiter-analogs. Therefore, we checked what fraction of our satellites could be detected with current and near-future instruments. Assuming the planet is already a transiting one, a sensitivity of 10 −5 would allow to detect 14% of moon-transits. To detect half of our population, a sensitivity of 2 × 10 −6 would be necessary. Since the population is relatively small in mass, their transit-timing-variation reaches above 10 seconds only in 6.2% of the case, posing a challenge to observe these moons.
We checked the difference made by having a 3D N-body model with 1D disc treatment, versus a fully 1D model (disc + dynamics). We found that the migration rates in the 3D case are lower, and hence the moons can grow to larger masses, and less frequently migrate all the way into the planet. Furthermore, in the fully 1D model, the scattering of satellites cannot be accounted for due to the limited dimensionality. In the 3D runs, there is a scattered, inner small moon population, that is missing from the fully 1D runs.
We also examined the difference between 1D and 2D treatment of the circumplanetary disc in the planet formation model. The dust accretion proceeds differently, since only the dust within the feeding zone can be accreted (and therefore depleted) in the 2D case, while in the 1D case the entire azimuthal ring is available for accretion. This leads to the fact that in the 2D case initially the moons are smaller and accreting with a lower rate in comparison to the 1D. As the moonlets grow, their feeding zone enlarges too, even in the 2D case. The differential velocity within the feeding zone will then provide more dust available for the moon to accrete in the 2D case, and by the end of the formation, the results on the satellite masses between the 1D and 2D simulations are small, never more than 15%. Figure A1. The figure shows the position of systems in a 1D projection got by a t-SNE approach as a function of the dust-to-gas ratio and the refilling timescale. Some correlation between this two parameters is easily visualized.

APPENDIX A: RANDOMIZED DEPENDENCE COEFFICIENT
In Section 3.3 an important parameter has been identified, i.e. Γ = ref / 2.2 , that is supposed to control the architecture of systems. Simulations with the same ratio between these two quantities will produce similar outcomes. The relation between them has been checked in Figure 11, 12, and 13, but it can also be formally found thanks to a procedure known as Randomized Dependence Coefficient (Lopez-Paz et al. 2013;Doran 2017). The RDC is a measure of nonlinear dependence between random variables based on the Hirschfeld-Gebelein-Renyi Maximum Correlation Coefficient (Hirschfeld 1935). This coefficient will be 1 if two variables are correlated, even if non-linearly, 0 if they are not.
Given that the initial number of seeds does not influence the results significantly, the idea is that therefore the outcomes of the model would depend on some combination of the two other parameters. Mathematically, some quantity that identifies the outcomes would be a function of this combination, i.e. where Q could be for example the total mass of satellites. Even better, if we refer to the framework presented in Section 3.3, this quantity could be the position of systems on one axis as computed by the t-SNE algorithm. In fact, even if t-SNE has been used to represent systems in a 2D space, it is also able to reduce the dimensions further and represent systems on a 1D line. That means, in practice, unfolding and projecting the tube-shape in Figure 13 on one dimension.
If this intuition is true and we choose the right value for , then the RDC between = ref and would be 1, or close to it. The procedure is then maximizing the RDC between and while varying , in order to find which combination of the two parameters is correlated the most with the results. It is already expected that a combination exists, given the results already shown in Section 3.3. In order to do that, we studied the problem in a log-space, given that the results spanned over many orders of magnitude. Then the correlation can be examined between = log 10 = log 10 + log 10 ref and the t-SNE 1D distance ( Figure A1 and A2). Figure A1 shows the position of systems as a function of the two parameters and ref , while Figure A2 shows the distance as a function of the for which the RDC has been maximized. In this case the algorithm has recognized a 95% non-linear correlation with = −0.45. This means that there is a non-linear correlation between the outcomes of the model and the parameter ∼ ref 2.2 .

APPENDIX B: OBSERVABILITY MODEL
In Section 3.5, we widely used the concept of detection probability for a given system around a Jupiter-like planet. First of all, anytime we considered this probability, we implied that that Jupiter-like planet is already detected and transiting. Otherwise the factor of ★ 7×10 7 5.2 ∼ 10 −4 has to be added in front of any probability (Perryman 2011).
Then, for each satellite in a given system we assigned a probability . Once we have calculated those, the probability of detecting at least one satellite in that system would be There are two effects that contribute to . The first one is trivial, i.e. we assign a probability of 1, if the transit depth of the satellite ( sat / ★ ) is bigger than the sensitivity, while we assign 0 in the other case. Still, if a satellite is in principle detectable, we have to decrease the probability because of geometrical factors. These factors actually represent the second effect, which has to do with the fact that, even if we detect a Jupiter-like planet, a satellite could be not visible because it is hiding behind or in front of the planet or because it is not transiting in front of the star (i.e. having high inclination).
In the first case, we refer to the scheme in Figure B1. There is a portion of the orbit of a satellite during which it is passing in front of Figure B1. Face-on configuration of the first geometrical effect. The satellite is transiting behind or in front of the planet during a portion of its orbit. Figure B2. Side-on scheme of the second geometrical effect. Even though the planet is transiting in front of the star, a satellite may spend part of its orbit out of the stellar disc.