The Role of Drag and Gravity on Dust Concentration in a Gravitationally Unstable Disc

We carry out three dimensional smoothed particle hydrodynamics simulations to study the role of gravitational and drag forces on the concentration of large dust grains (St>1) in the spiral arms of gravitationally unstable protoplanetary discs, and the resulting implications for planet formation. We find that both drag and gravity play an important role in the evolution of large dust grains. If we include both, grains that would otherwise be partially decoupled will become well coupled and trace the spirals. For the dust grains most influenced by drag (with Stokes numbers near unity), the dust disc quickly becomes gravitationally unstable and rapidly forms clumps with masses between 0.15 - 6 Earth masses. A large fraction of clumps are below the threshold where runaway gas accretion can occur. However, if dust self-gravity is neglected, the dust is unable to form clumps, despite still becoming trapped in the gas spirals. When large dust grains are unable to feel either gas gravity or drag, the dust is unable to trace the gas spirals. Hence, full physics is needed to properly simulate dust in gravitationally unstable discs. Dust trapping of large grains in spiral arms of discs stable to gas fragmentation could explain planet formation in very young discs by a population of planetesimals formed due to the combined roles of drag and gravity in the earliest stages of a disc's evolution. Furthermore, it highlights that gravitationally unstable discs are not just important for forming gas giants quickly, it can also rapidly form Earth mass bodies.


INTRODUCTION
In recent years, the Atacama Large Millimeter/submillimeter Array (ALMA) has revealed that many of the brightest, and presumably the most massive, protoplanetary discs are highly structured.The most common of these sub-structures are rings & gaps (ALMA Partnership et al. 2015;Andrews et al. 2016Andrews et al. , 2018;;Fedele et al. 2018;Huang et al. 2018a;Dipierro et al. 2018;Long et al. 2018;Booth & Ilee 2020).A natural explanation for the origin of ring & gap structures is through planet-disc interaction (see reviews by Kley & Nelson 2012;Paardekooper et al. 2022).The evidence of planets being the cause has been strengthened through studies of the gas kinematics (Perez et al. 2015;Pinte et al. 2018Pinte et al. , 2019Pinte et al. , 2020;;Calcino et al. 2022;Pinte et al. 2023a).
The discs with observed ring & gap structure are young, with some being less than 1Myr old (ALMA Partnership et al. 2015;Fedele et al. 2018;Dipierro et al. 2018;Sheehan & Eisner 2018; Segura-Cox et al. ★ E-mail: sahl.rowther@leicester.ac.uk 2020).Forming massive planets at these locations is challenging at such young ages.Many of the potential planets inferred from the rings & gaps and kinematics are massive (0.1 − 10  Jup ) and far from the central star at tens or even hundreds of au (Lodato et al. 2019;Pinte et al. 2023b).Traditional planet formation models such as core accretion struggle to form massive planets at wide orbital separations so quickly.Hence, a natural question to ask is how these putative planets could have formed so quickly.
An alternate planet formation mechanism is through gravitational instability.In the earliest stages of a disc's evolution, it is expected to be massive enough to develop gravitational instabilities in the form of spiral structures.Although this phase is short-lived, there is observational evidence of such discs (Pérez et al. 2016;Huang et al. 2018b;Paneque-Carreño et al. 2021), which is also supported by theoretical predictions (Meru et al. 2017;Hall et al. 2020;Longarini et al. 2021).In this scenario, if the disc is massive enough, it can fragment and form giant planets (Boss 1997) on dynamical timescales (Gammie 2001).An issue is that although the initial mass of these fragments is in the planetary regime, they quickly accrete a lot of mass and become brown dwarfs ( >13 Jup ) as they migrate (Stamatellos & Inutsuka 2018).However, there is another path to forming planets if the disc is stable to fragmentation, but still massive enough to form spirals.The spirals are regions of high pressure where dust can become trapped and grow to form planetesimals which become the seeds of planet formation (Rice et al. 2004(Rice et al. , 2006;;Gibbons et al. 2012;Booth & Clarke 2016;Elbakyan et al. 2020;Baehr & Zhu 2021;Baehr et al. 2022;Longarini et al. 2023a,b).On the other hand, it has also been shown that gravitational instability can inhibit the concentration of dust (Walmswell et al. 2013;Riols et al. 2020).The differing conclusions could be due to the widely different methodologies which range from 2D local shearing boxes to 3D global simulations, and from using dust particles in the test particle regime to feeling gravitational acceleration from both gas and itself.Thus, the aims of this study (using 3D global simulations) are to determine what drives the dynamics of large dust grains in these spirals; the drag force, or gravity?
In this study we perform three-dimensional global numerical simulations to investigate the role of drag and gravity on the dust dynamics in a gravitationally unstable protoplanetary disc.In §2 we describe the simulations presented in this study.In §3 we present the results of dust concentration in the spiral arms, and roles of drag and gravity.The limitations and implications of our work are discussed in §4.

MODEL
We use Phantom, a smoothed particle hydrodynamics (SPH) code developed by Price et al. (2018) to perform the suite of simulations presented here.SPH codes are favourable for scenarios with high density contrasts.Thus, Phantom is well suited for the problem at hand and has been previously used for dusty simulations (Laibe & Price 2012a,b;Dipierro et al. 2015;Price & Laibe 2020) and for simulations involving self-gravity (Rowther & Meru 2020;Cadman et al. 2020;Rowther et al. 2020Rowther et al. , 2022Rowther et al. , 2023)).For the simulations in this work, the dust and gas are modelled as two separate fluids (Laibe & Price 2012a,b).

Disc setup
All discs are modelled using 2 × 10 6 gas particles and 2.5 × 10 5 dust particles between  in = 4au and  out = 100au where the fiducial disc has a mass of 0.2 ⊙ disc around a 1 ⊙ star.A sink particle (Bate et al. 1995) is used to model the central star.The accretion radius of the central star is set to be equal to the disc's inner boundary,  in .Particles within 0.8 times the accretion radius are accreted immediately without any checks.The surface density profile Σ is given by where Σ 0 = 1.1×10 3 g cm −2 , and   = 1 − √︁  in / is the factor used to smooth the surface density at the inner boundary of the disc.The initial temperature profile is expressed as a power law where  0 is set such that the disc aspect ratio / = 0.05 at  =  0 =  in .The energy equation for the gas particles is where we assume an adiabatic equation of state with  = 5/3, and  is the specific internal energy,  is the pressure,  is the density and  is the velocity.The subscript 'g' refers to gas particles.We do not include drag heating from back-reaction of the dust on the gas.The first term on the RHS is the d work, and Λ shock is a heating term that is due to the artificial viscosity used to correctly deal with shock fronts.The final term controls the cooling rate in the disc.Here the cooling time is straightforwardly implemented to be proportional to the dynamical time by a factor of , where Ω is the orbital frequency, and  = 15.The Cullen & Dehnen (2010) switch detects any shocks that form and generates the correct artificial viscosity response.The linear artificial viscosity parameter  AV varies depending on the proximity to a shock.Close to the shock, it takes a maximum of  max = 1, and a minimum of  min = 0.1 far away.The quadratic artificial viscosity coefficient  AV is set to 2 (see Price et al. 2018;Nealon et al. 2015).

Dust dynamics
The dust dynamics is governed by gravitational and drag forces from interactions with the gas.The acceleration of the dust particles is given by where the subscript 'd' refers to dust particles.The first term on the right hand side is the drag force which acts to eliminate any velocity difference between the dust and gas.The second term on the right hand side of Equation 6 is the gravitational acceleration due to interactions between dust and gas particles; the third term is the dust self-gravity; and the final term is the gravitational acceleration from the central star.
An important timescale for the drag force is the stopping time,  s , the timescale at which the differential velocity between the dust and gas decays.If the stopping time is short, the dust velocities quickly decay to the gas velocities resulting in the dust becoming coupled to the gas.If it's long, the dust velocities will remain different to the gas velocities, and hence the dust structure will be different from the gas structure.The stopping time is related to the drag coefficient  (Laibe & Price 2012b), by where  =  g +  d .However, to decrease computational expense when dust clumps form, we instead approximate  to be equal to  g only when determining the minimum timestep required for evolving the particle.The effectiveness of the drag force can be described by the Stokes number St, which is defined by the ratio of the stopping time to the orbital timescale, The simulations in this study fall in the Epstein drag regime where the stopping time is approximated by where  = 5/3 is the adiabatic index,  s is the sound speed,  grain is the grain size,  grain = 3g/cm 3 is the intrinsic grain density, and  is a correction for supersonic drift velocities given by Kwok (1975) as where Δ  ≡ |  −   |.For small Stokes numbers (< 1), the stopping time is sub-orbital, and the drag force quickly eliminates any velocity difference between the dust and gas.However, for large Stokes numbers (> 1), it can take multiple orbits before the velocity difference decays.Hence for a dynamic disc where the substructure continuously evolves -as is the case for a gravitationally unstable disc -it becomes difficult for the dust velocities to decay to the gas velocities.

Defining Toomre Q for dust
The Toomre Q parameter (Toomre 1964) gives a measure of how gravitationally unstable a disc is, and is defined as where  is the gravitational constant,  s is the sound speed, and Σ is the surface density.A disc becomes more gravitationally unstable as the surface density increases, or as the sound speed decreases.Although the above equation is generally used for the gas, we calculate the Toomre Q of dust using the effective sound speed and surface density of the dust.The effective sound speed of the dust particles  dust s can be obtained by calculating their velocity dispersion as where the calculation is summed over all dust particles that are neighbours of particle .The mass and density are given by   and  d,  , and  d, and  d,  are the magnitudes of the velocity of dust particles  and , respectively.The smoothing kernel and length are given by    and ℎ  , respectively.

The suite of simulations
Initially, the disc is modelled as a gas only simulation until spiral structures have developed.After 4 outer orbits, the dust is added uniformly as a separate dust disc as if it is  = 0 where no spirals are present (see bottom left panel of Figure 1).In the fiducial simulation, 50 cm sized dust grains are used.This corresponds to an initial average Stokes number of 4 at  = 50au.The effectiveness of the drag force is dependent on  s , and thus the size of the dust grains, where the larger the dust grain, the less efficient the drag force becomes.Since the drag force acts to eliminate any velocity difference between the dust and gas, a weakening drag force leads to decoupling of the dust.Hence, to determine when dust should stop tracing the gas spirals in a gravitationally unstable disc, three additional simulations were modelled with dust sizes of 150, 500, and 5000 cm.This corresponds to initial average Stokes numbers of 12, 40, and 400, respectively.
To investigate the importance of gravitational and drag forces, the simulation with 50 cm sized grains at  = 4 outer orbits is used as the initial condition for four further simulations.The initial conditions can be seen in the left column of Figure 1.Except for the reference simulation, each simulation turns off certain physics to see how the dust evolves in its absence.The four simulations are Full Physics -The reference simulation.The evolution of the fiducial disc is continued as normal where the dust feels both the gravity of gas and dust, as well as the drag force.
No Dust Self-Gravity -The dust no longer feels the gravity of other dust particles.It only feels the gas gravity and the drag force.Computationally, dust self-gravity is turned off by decreasing the dust particle mass by eight orders of magnitude.This makes the third term in Eq 6,  dust−dust negligible but does not affect the other terms.
Drag Only -The dust is no longer affected by gravity.It only feels the drag force.The second term in Eq 6,  dust−gas is set to 0, in addition to the decreased dust particle mass.However, since the gas still feels its own self-gravity, the velocity of the gas particles is increased due to the gas disc itself.If the dust has no knowledge of the gas disc at all, the dust particles will orbit at a slower velocity relative to the gas, resulting in outward migration.Thus, the dust velocity includes a correction for the enclosed mass in the Keplerian velocity.The dust therefore has no local knowledge of any gravitational potential, but only the global gas disc.

Gravity Only -
The drag force has no effect on the dust.Only gas and dust gravity acts on the dust.Here, only the first term (the drag force) in Eq 6 is set to 0. This is a similar scenario to the simulations in Walmswell et al. (2013).
In all simulations, the dust always feels the gravitational acceleration from the central star (the final term in Eq 6), and the gas always feels its own self-gravity.

Clump finding
We use the clump finding algorithm introduced in Wurster & Bonnell (2023).In a clump, the densest particle is considered to be 'lead'.All other particles in the clump are considered to be 'members'.Briefly, we first sort dust particles by density.Then, the densest particle with  >  lead = 2 × 10 −13 g cm −3 becomes the lead member of the first clump.For successively less dense particles, we determine if their smoothing length overlaps with an existing clump and if the particle is bound to the clump.If so and if  >  member = 10 −13 g cm −3 , then it is added to the clump and the clump properties are updated to account for the new member.If not and  >  lead , then it becomes the lead member of a new clump.When a particle's smoothing length overlaps with two clumps, or if two clumps overlap, we determine if they are bound, and if so, we merge the clumps.These values are chosen to help the clump-finding algorithm distinguish clumps from the spiral arm they reside in.The impact of our choices do not meaningfully affect the analysis of the clumps.
Since we are evaluating clump membership by inspecting particles in order of decreasing density rather than spatial proximity to a clump, we repeat this process iteratively until the number of clumps and their properties have converged.This is required since each iteration will find more distant particles whose smoothing length will overlap with the clump.Only dust particles have been considered for the clump finding algorithm.Clumps with less than 58 particles (the average number of neighbours) are considered to be unresolved, and are removed.

Dust concentration and clumping in the spiral arms
Figure 1 shows the density evolution of gas (top) and 50cm sized dust grains (bottom) in a gravitationally unstable 0.2 ⊙ disc at three moments in time.The left column shows the gas and dust at a time when the dust is initialised.The gas disc is stable to fragmentation.The middle and right columns show the simulation 0.44 and 0.88 outer orbits later, respectively.The 50cm dust grains initially have an average Stokes number, St ∼ 4. In this regime, dust particles are expected to drift efficiently into pressure maxima (Weidenschilling 1977).The spiral arms due to gravitational instability are regions of pressure maxima.They are also regions of gravitational potential minima which also attracts the dust.Hence, dust rapidly drifts towards and concentrates in the spiral arms in the middle and right panels of Figure 1.The dust spirals are narrower compared to the gas spirals.
Figure 2 shows the dust-to-gas mass ratio (top), the sound speed of the dust relative to the gas (middle), and Toomre Q of the dust disc (bottom).Given that dust is only 1% of the entire disc mass, it is perhaps surprising that there is enough dust mass to become gravitationally unstable.Especially since the initial  dust min ∼ 90, which is far from the conditions required for gravitational instability to occur ( < 1.7, Durisen et al. 2007).As the disc has evolved, the local dust-to-gas mass ratio has increased from its initial value of 0.01.In spiral arms the dust-to-gas mass ratio has increased by an order of magnitude over the initial value.The sound speed of the dust is also lower relative to the gas by nearly an order of magnitude.Both of these factors explain why the dust disc collapses into clumps despite the initial dust profile being safely in the gravitationally stable regime.The increased surface density and lower sound speed relative to the gas both contribute to the dust disc becoming locally gravitationally unstable in the spiral arms, with  dust < 1.7 as shown in the bottom panel of Fig 2.

Radial Drift
In contrast to the gas, the dust is a pressureless fluid which results in their velocities being different.In a smooth axisymmetric disc the drag force exerts a torque on the dust particles, resulting in dust migrating inwards towards the central star.However, in these simulations the gas structure is not axisymmetric.The spiral arms due to gravitational instabilities are local pressure maxima, and dust therefore drifts towards them.This can be seen in Figure 3, which shows the cross section slice of the radial velocity,   , of the dust in the  = 0 plane.After the initial drift towards the spiral arms, the dust motion is along the spirals.

Larger dust sizes
Even the smallest grain size (50cm) in this study has Stokes numbers above unity.In non-self-gravitating discs in this regime, the dust is partially decoupled from the gas.However, as seen from Figure 1, in a gravitationally unstable disc, the dust still traces the gas structure.The dust is not perfectly coupled as evident from the enhanced dustto-gas mass ratios.
Figure 4 shows the cross sectional slices of the dust density for all four dust sizes.Spiral structures are visible in all four panels, although to different extents.As the dust size increases, the dust  The structure of the 150 cm sized grains is the most similar to the gas structure.This can be observed by the peak of the PDF at  dust s ∼  gas s .Essentially, the dust and gas act similarly since the similar sound speeds cause the response of both fluids to any perturbations to be similar.Thus, the dust structure is closely coupled to the gas.At larger grain sizes, the drag force becomes less effective in damping velocity variations resulting in  dust s >

The role of self-gravity and drag
It's apparent that drag plays an important role in shaping the dust as seen by the evolution of different sized dust grains.The role of gas gravity and dust self-gravity remains less clear.Here we focus on the 50 cm sized grains, where the drag force is most efficient for this simulation, making it easier to study the role of self-gravity compared to the drag force.
Figure 6 shows the cross sectional slices of the dust density 0.88 outer orbits after the simulation is resumed with different combinations of physics turned off.The reference simulation continues  to evolve as expected.The dust remains trapped in the spiral arms, with numerous clumps.When the dust cannot feel the gravity of other dust particles (the "No Dust Self-Gravity" simulation), clumps cannot form.The dust evolves similarly to the reference simulation otherwise with dust still being trapped by the gas spirals.When dust cannot feel the gas gravity either (the "Drag Only" simulation), the continued presence of the drag force keeps the dust in thin filaments.However, as the dust no longer feels the gravitational forces due to the gas, it no longer traces the spiral structure (also see Figure A2).Drag alone isn't sufficient to trap dust and form clumps. Finally, Figure 6.Cross section slices of dust density in the  = 0 plane 0.88 outer orbit after physics has modified.The four simulations are: the reference simulation with physics as normal (Full Physics -top left panel), the dust no longer feels the gravitational force of other dust particles (No Dust Self-Gravity -top right panel), the dust no longer feels gravitational acceleration from either the gas or dust (Drag Only -bottom left panel), and the dust no longer feels the drag force (Gravity Only -bottom right panel).The dust densities are lower in the top right and bottom left panels because the dust particle mass is decreased by eight orders of magnitude to mimic turning off dust self-gravity.The loss of dust self-gravity results in clumps being unable to form despite the dust still becoming trapped in the gas spirals.Neither the drag force nor gravity is able to couple the dust to the gas.Both are necessary for the dust to trace the gas spirals.See the top right panel in Figure 1 for the gas distribution.The spiral structures visible in the bottom panels eventually disappear in the absence of drag/gravity (see Figure A2).
when the dust cannot feel the drag force, the spiral structures visible initially form due to gravitational interactions between dust and gas.However, without drag, the spirals are not maintained and disappear after a couple of outer orbits.The dust behaves similarly to the largest dust size in §3.3; the dust cannot form clumps and does not trace the gas spirals.
To illustrate the role gravity and drag have on the dust dynamics, consider the plot of the PDF of  dust s / gas s in Figure 7 along with the two timescales described in §2.2: the stopping time; and the orbital time.For the dust sizes (50cm, St ≈ 4) in Figure 6, the stopping time is larger than the orbital time.The spiral structures are constantly moving at the orbital velocity.The longer timescale of the drag force causes the dust to decouple when it only feels the drag force (bottom left panel in Figure 6).The dust velocities are unable to decay to the gas velocities within the time they cross the gas spirals.However, the dust spirals are also regions of gravitational potential minima.If the dust only feels gravitational forces (bottom right panel in Figure 6), the dust initially forms spiral structures but quickly decouples from the gas as there is no longer any force eliminating any differential velocity between the dust and gas resulting in  dust s >  gas s .This is essentially what is seen with larger dust sizes in Section 2.2.While neither drag nor gravity on their own is enough to trap dust, their combined influence enables the dust to keep up with the gas spirals (top right panel in Figure 6).
However, this analysis is specific to the large grain sizes studied here (St > 1).Smaller grain sizes with St < 0.1 will still be able to trace the spirals due to the drag force alone: the short stopping time results in the dust velocities decaying to the gas velocities before the gas spirals have moved significantly (Baehr & Zhu 2021).
As mentioned earlier, all dust sizes here are in the regime where dust is expected to be decoupled from the gas (St > 1).So it's perhaps a bit surprising that the dust was still able to trace the gas spirals.The gravity of the gas is usually negligible in most studies of dusty protoplanetary discs which are gravitationally stable.These results show that in a gravitationally unstable disc, the combined effect of gas gravity and drag is necessary for large dust grains to be well coupled , for all four simulations with varying physics as described in Figure 6 and in §2.4.When the dust only feels the drag force or gravity, it results in an increase in velocity dispersion resulting in dust decoupling from the gas.
to the gas and trace the spirals.For the dust to form clumps, dust self-gravity is necessary.Although including full physics results in a simulation that is at least an order of magnitude slower to compute (without drag), it is important to do so, particularly for studying the formation of dust clumps.

Clump formation
The top panel of Figure 8 shows a histogram of all resolved clump masses at  = 4.88 outer orbits for the reference (Full Physics) simulation where the dust feels the gravity of both the gas and dust, and the drag force.Figure A1 shows the locations of all resolved clumps.The clumps range from 0.15 ⊕ (the resolution limit) to ∼6 ⊕ , with a high majority of the clumps being less than the critical mass (∼5 ⊕ ) to undergo runaway gas accretion (Papaloizou & Terquem 1999;Rice & Armitage 2003).The critical mass is dependent on the disc parameters, and solid accretion rate and composition (see review by Drążkowska et al. 2023).The lower panel of Figure 8 shows the total mass found in gravitationally bound clumps as a function of time Clump formation is rapid with roughly a Jupiter mass of dust accumulating in clumps within half an outer orbit.The total mass in clumps steadies after this initial rapid build-up of clumps once roughly half the dust is in clumps.

Limitations
The main limitations of this study involve how the dust is modelled.The dust used in the simulations in this study is a single fixed size.While this simplifies the simulations, it is not representative of reality where dust is made up of a mixture of different sizes and can grow.Assuming all the dust is of the same size (and therefore all the dust mass is in a single grain size, rather than a power-law distribution) is the optimal scenario for the dust disc to become gravitationally unstable and form clumps. various dust sizes, then there might not be enough mass for dust of a given size to become gravitationally unstable.With a range of dust sizes, and hence Stokes numbers, the effectiveness of dust trapping will be varied.Not all dust will be optimally trapped, thus resulting in reduced dust density enhancements and increasing .However, including dust growth could aid dust clumping.Small dust grains are well coupled to the gas and collect in the gas spirals (Baehr & Zhu 2021;Gibbons et al. 2012).In these regions, Elbakyan et al. (2020) showed that dust can rapidly grow from micron sized grains to a few centimetres with St ∼ 0.3.If enough dust can grow to this size where they are most strongly influenced by the gas spirals, then the dust might still be able to become gravitationally unstable and form clumps.
Our simulations have not included dust back-reaction onto the gas.With the enhancement of the dust-to-gas mass ratios as seen in Figure 2, this could become relevant.However, based on shearing box simulations in Gibbons et al. (2014) and Baehr et al. (2022), the qualitative nature of the role of drag and gravitational forces as described in this study is likely to be unchanged if back-reaction was included.However, the mass of the clumps measured would be impacted.Baehr et al. (2022) found clump masses to be smaller, but more numerous when back-reaction was included.
These simulations also made use of a simplifying approximation to the stopping time.This simplification allowed for the drag force term in Eq 6 to be independent of the dust density which enables the simulations to evolve for a longer time after dust clumps are formed.Our results of clump formation still compare well with Rice et al. (2006) which is similarly as rapid.While their simulations stopped when 15% of dust was in clumps, the amount of dust in clumps was still rising.As we could evolve further, we find that the dust continued to accumulate in clumps until ∼50% of the dust was in clumps.We find the clump masses (∼0.1 − 5 ⊕ ) are slightly higher than the results in Baehr et al. (2022), but similar to Longarini et al. (2023a,b), although this could be due to the lower Stokes numbers in the former's simulations.Additionally they utilised a fixed stopping time, whereas we used a fixed particle size which results in the stopping time varying with the particle location.A fixed particle size has been shown to result in stronger dust concentrations (Shi et al. 2016).

Planet formation
The simulations here show that planet formation in gravitationally unstable discs is not restricted to gas giants.The spiral arms of gravitationally unstable discs provide conditions that favour dust concentration which rapidly collapses to form numerous Earth mass clumps.The short timescale at which the Earth mass clumps form could help explain how some of these discs exhibit ring & gap structure (often associated with evidence of planet formation) despite being less than a million years old (ALMA Partnership et al. 2015;Sheehan & Eisner 2018;Segura-Cox et al. 2020).The clumps formed in gravitationally unstable discs could become the seeds of the planets that go on to carve rings & gaps in ALMA observations.This work adds to the idea that gravitationally unstable discs may form planets lower in mass than giant planets (Deng et al. 2021).
The simulations in this study highlight that dust trapping in the spiral arms of non-fragmenting gravitationally unstable discs is a viable mechanism for forming planetesimals and/or planetary cores.Most discs are expected to go through a gravitationally unstable phase during their evolution.Thus, planet formation in very young discs could be explained by a population of planetesimals formed in the earliest stages of a disc's evolution (Nixon et al. 2018).However, the runtime of these simulations is restricted as the minimum timestep required for accurate evolution decreases as dust particles get closer together as it clumps.Future global simulations with dust growth are necessary to evolve these simulations further to investigate whether the clumps formed here remain as Earth to Super-Earths or whether they eventually become giant planets.Recently, Baehr (2023) have explored this question, though not with global simulations.

CONCLUSION
We perform 3D SPH simulations to investigate the role of gravitational and drag forces on dust concentration in gravitationally unstable protoplanetary discs.In the regime studied in this work (large dust grains with Stokes numbers greater than unity), our simulations show that both drag and gravity are necessary for dust to be well coupled to the gas and trace the spirals.Neither on their own is capable.
If the dust is strongly influenced by the drag force, the density increases while the sound speed of the dust decreases.The combined effect results in the dust becoming gravitationally unstable and forming gravitationally bound clumps.The speed at which the dust forms Earth to Super-Earth mass clumps is rapid.Within half an outer orbit, roughly half the dust disc is in the form of clumps.After which, clump formation ceases (as the disc runs out of material to form more).While dust is still trapped in the gas spirals without dust self-gravity, it is unable to collapse to form clumps. Neglecting either gas gravity or drag results in the dust being unable to trace the gas spirals.
The results presented here show that planet formation through gravitational instability may not be limited to just giant planets.If the clumps formed remain under the critical mass for runaway gas accretion (∼5 ⊕ ), then gravitational instability in discs stable to gas fragmentation could form Earth to Super-Earths very quickly through efficient dust trapping in the gas spirals.

Figure 1 .
Figure 1.Cross section slices of density in the  = 0 plane of gas, i.e. the midplane, (top) and 50cm sized dust grains (bottom) at  = 4, 4.44, and 4.88 outer orbits (from left to right).The dust is initialised at  = 4 orbits after spiral structures have formed in the gas.Dust rapidly concentrates in the spirals becoming gravitationally unstable and forming clumps as seen at  = 4.44, and 4.88 outer orbits.At  = 4.88 outer orbits, a few of the dust clumps have become massive enough to cause spiral wakes which are visible in both the gas and dust.

Figure 2 .
Figure 2. The top and middle panels show the cross section slices of the dustto-gas mass ratio and the sound speed of dust relative to the gas, respectively in the  = 0 plane.The bottom panel shows Toomre  of the dust disc.The combined effect of the increased dust enhancement in the spirals and lower sound speed results in the dust disc becoming locally unstable in the spiral arms (red regions).

Figure 3 .
Figure 3. Cross section slice of the radial velocity   of the dust in the  = 0 plane.The dust drifts towards the gas spirals which results in dust becoming trapped.The presence of gravitational instabilities could slow down radial drift of dust towards the central star.The orange and purple regions correspond to outward and inward dust motion, respectively.We can clearly see drift towards the spiral arms from both directions.

Figure 5
shows the probability density function (PDF) of the sound speed of dust relative to the gas,  dust s / gas s , for all four dust sizes.To compare  dust s and  gas s with each other, they were first individually interpolated onto a grid.The  dust s increases as grain size increases.Both the more efficient dust trapping and decreased dust sound speed are why only the smallest grain sizes are prone to clumping.
gas s .Hence the dust structure becomes more uncorrelated with the gas with increasing dust size.For the 5000 cm sized grains, the spiral structures visible in the bottom right panel of Fig 4 disappear after a few outer orbits.The PDF of the 50 cm sized grains peaks at  dust s visible which allows dust to collapse further into clumps.

Figure 4 .
Figure 4. Cross section slice of the dust density (left panels) and the dust-to-gas mass ratio (right panels) in the  = 0 plane 1.5 outer obits after the dust was added.Each panel represents a different size: 50 cm (top left), 150 cm (top right), 500 cm (bottom left), and 5000 cm (bottom right).The corresponding Stokes numbers are roughly 4, 12, 40, and 400 respectively.As the dust size increases, dust trapping becomes less efficient due to the decreasing influence of the drag force.Only the smallest two dust sizes become gravitationally unstable and form clumps.

Figure 5 .
Figure 5.A probability density function (PDF) of the sound speed of dust relative to the gas,  dust s / gas s (left panel), for all four dust sizes as described in Figure4.The drag force becomes less effective with increasing dust size.Hence, the dust motion becomes uncorrelated with the gas as seen by the increasing velocity dispersion with increasing dust size.The right panels show the Toomre Q of the dust disc for all four dust sizes.The less efficient dust trapping with increasing dust size results in the dust becoming more stable.

Figure 7 .
Figure7.A probability density function (PDF) of the sound speed of dust relative to the gas,  dust s / gas s , for all four simulations with varying physics as described in Figure6and in §2.4.When the dust only feels the drag force or gravity, it results in an increase in velocity dispersion resulting in dust decoupling from the gas.

Figure 8 .
Figure 8.The top panel shows a histogram of the clump mass at  = 4.88 outer orbits.A majority of the clumps are less than a few Earth masses, and hence might not be massive enough to become giant planets through runaway gas accretion.The bottom panel shows the time evolution of the total mass found in clumps.Within half an outer orbit 60% of the dust disc (400 ⊕ ) is found in clumps as shown by the black line (Full Physics simulation).Only resolved clumps are plotted.

Figure A2 .
Figure A2.Cross section slices of dust density in the  = 0 plane showing the evolution of the bottom two panels in Figure 6.The two simulations are where the dust no longer feels gravitational acceleration from either the gas or dust (Drag Only -left panel), and the dust no longer feels the drag force (Gravity Only -right panel).The spirals that were present in Fig 6 eventually disappear if the dust only feels drag or gravity.The combined influence of drag and gravity is required for dust to remain trapped in spirals.