  The origins of the Galilean satellites - namely Io, Europa, Ganymede, and
Callisto - is not fully understood yet. Here we use N-body numerical
simulations to study the formation of Galilean satellites in a gaseous
circumplanetary disk around Jupiter. Our model includes the effects of pebble
accretion, gas-driven migration, and gas tidal damping and drag.
Satellitesimals in our simulations first grow via pebble accretion and start to
migrate inwards. When they reach the trap at the disk inner edge, scattering
events and collisions take place promoting additional growth. Growing
satellites eventually reach a multi-resonant configuration anchored at the disk
inner edge. Our best match to the masses of the Galilean satellites is produced
in simulations where the integrated pebble flux is 1e-3 MJ. These simulations
typically produce between 3 and 5 satellites. In our best analogues, adjacent
satellite pairs are all locked in 2:1 mean motion resonances. However, they
have also moderately eccentric orbits (0.1), unlike the current real
satellites. We propose that the Galilean satellites system is a primordial
resonant chain, similar to exoplanet systems as TRAPPIST-1, Kepler-223, and
TOI-178. Callisto was probably in resonance with Ganymede in the past but left
this configuration - without breaking the Laplacian resonance - via divergent
migration due to tidal planet-satellite interactions. These same effects
further damped the orbital eccentricities of these satellites down to their
current values (0.001). Our results support the hypothesis that Io and Europa
were born with water-ice rich compositions and lost all/most of their water
afterwards. Firmer constraints on the primordial compositions of the Galilean
satellites are crucial to distinguish formation models.
