The Limited Role of the Streaming Instability during Moon and Exomoon Formation

It is generally accepted that the Moon accreted from the disk formed by an impact between the proto-Earth and impactor, but its details are highly debated. Some models suggest that a Mars-sized impactor formed a silicate melt-rich (vapor-poor) disk around Earth, whereas other models suggest that a highly energetic impact produced a silicate vapor-rich disk. Such a vapor-rich disk, however, may not be suitable for the Moon formation, because moonlets, building blocks of the Moon, of 100 m–100 km in radius may experience strong gas drag and fall onto Earth on a short timescale, failing to grow further. This problem may be avoided if large moonlets (≫100 km) form very quickly by streaming instability, which is a process to concentrate particles enough to cause gravitational collapse and rapid formation of planetesimals or moonlets. Here, we investigate the effect of the streaming instability in the Moon-forming disk for the first time and find that this instability can quickly form ∼100 km-sized moonlets. However, these moonlets are not large enough to avoid strong drag, and they still fall onto Earth quickly. This suggests that the vapor-rich disks may not form the large Moon, and therefore the models that produce vapor-poor disks are supported. This result is applicable to general impact-induced moon-forming disks, supporting the previous suggestion that small planets (<1.6 R ⊕) are good candidates to host large moons because their impact-induced disks would likely be vapor-poor. We find a limited role of streaming instability in satellite formation in an impact-induced disk, whereas it plays a key role during planet formation.


INTRODUCTION 1.Review of the Moon-formation hypothesis
The origin of the Earth's Moon is a long-standing open problem in planetary science.While it is accepted that the Moon formed from a partially vaporized disk generated by a collision between the proto-Earth and an impactor approximately 4.5 billion years ago (the "giant impact hypothesis") (Hartmann & Davis 1975;Cameron & Ward 1976), the details, such as the impactor radius and velocity, are actively debated.According to the canonical model, the proto-Earth was hit by a Marssized impactor (Canup & Asphaug 2001).This model can successfully explain several observations, such as the mass of the Moon, the angular momentum of the Earth-Moon system, and the small lunar core.Moreover, this model could potentially explain the observation that the Moon is depleted volatiles, if they escaped during the impact or subsequent processes (see reviews by Halliday 2004;Canup et al. 2023).However, this model fails to explain the Earth and Moon's nearly identical isotopic ratios (e.g., Si, Ti, W, O, Armytage et al. 2012;Zhang et al. 2012;Touboul et al. 2015;Thiemens et al. 2019Thiemens et al. , 2021;;Kruijer et al. 2021).This is because the disk generated by an oblique Mars-sized impactor is primarily made of the impactor materials (Canup 2004;Kegerreis et al. 2022;Hull et al. 2024), which likely had different isotopic ratios from Earth, unless the inner solar system was well mixed and homogenized (Dauphas 2017) or equilibrated in the disk phase (Pahlevan & Stevenson 2007;Lock et al. 2018).In contrast, more energetic impact models may solve this problem by mixing the proto-Earth and impactor.These energetic models include an impact between two half Earth sized objects (Canup 2012), an impact between a rapidly rotating proto-Earth and a small impactor (Cúk & Stewart 2012), and general impact events that involve high energy and high angular momentum, which creates a doughnut-shaped vapor disk connected to the Earth (a so-called Synestia) (Lock et al. 2018).In these scenarios, the disk and the proto-Earth compositions could have been mixed well, potentially solving the isotopic problem.Other proposed scenarios include (a) the Moon formation by multiple small impacts (Rufu et al. 2017), (b) hit-andrun collisions that mix the proto-Earth and the Moon more than the canonical model does (Asphaug et al. 2021), and (c) an impact between the proto-Earth with a surface magma ocean, which leads to more contribution from Earth to the disk than the canonical scenario (Hosono et al. 2019).
All of these models have potential explanations for the isotopic similarity, but each model faces challenges to explain other constraints.For example, the half-Earths, high angular momentum, and high energy models predict a much higher angular momentum of the Earth-Moon system than that of today (Cúk & Stewart 2012;Canup 2012).The excess of the angular momentum may or may not be removed easily (Cúk & Stewart 2012;Cúk et al. 2016;Ward et al. 2020;Rufu et al. 2017;Rufu & Canup 2020;Ćuk et al. 2021).Another constraint may come from the Earth's interior; the deep portion of the Earth still retains solar neon (Yokochi & Marty 2004;Williams & Mukhopadhyay 2019) and helium isotopic ratios (Williams et al. 2019), which are distinct from the rest of the mantle.These noble gasses may have been delivered to Earth's mantle during the planet formation phase.This suggests that Earth's mantle developed chemical heterogeneity during Earth's accretion, when the protoplanetary disk was still around, and that the Earth's mantle has never been completely mixed, even by the Moon-forming impact.Previous work suggests that preservation of the mantle heterogeneity can be explained by the canonical model, but not by the energetic models, because these energetic impacts tend to mix the mantle (Nakajima & Stevenson 2015).It should be noted, however, the possibility that the neon and helium heterogeneity developed after the Moon-formation has been discussed (Bouhifd et al. 2020), assuming delivery of these noble gases to the core was efficient.
Moreover, due to recent analytical capability, small isotopic differences between Earth and the Moon have been observed (e.g., K, O, W, V, Cr, and others, Wang & Jacobsen 2016;Wiechert et al. 2001;Young et al. 2016;Touboul et al. 2015;Kruijer et al. 2015;Thiemens et al. 2019;Nielsen et al. 2021;Sossi et al. 2018).Some isotopic difference, such as the Moon's enrichment of heavy K isotopes compared to those of Earth, can be explained by isotopic fractionation during the Moon accretion process due to liquid-vapor phase separation (Wang & Jacobsen 2016;Nie & Dauphas 2019;Charnoz et al. 2021).Observed small oxygen isotopic differences between Earth and the Moon may suggest that the Moon still records the impactor component (Cano et al. 2020).This suggestion may not be compatible with the energetic impact models because these impacts are so energetic that Earth and the Moon would be efficiently mixed.Whether the proposed models for the lunar origin can explain the observation that the Moon is depleted in volatiles is actively debated (Canup et al. 2015;Dauphas et al. 2015;Nakajima & Stevenson 2018;Dauphas et al. 2022;Nie & Dauphas 2019;Sossi et al. 2018;Charnoz et al. 2021;Halliday & Canup 2022;Canup et al. 2023).

Gas drag problem in a vapor-rich disk
Another key constraint that has not been discussed until recently is the vapor mass fraction of the disk (VMF), which sensitively depends on the impact model.Less energetic impacts, such as the canonical model and the multiple impact model, produce relatively small vapor mass fractions (VMF∼ 0.2 for the canonical model, Nakajima & Stevenson 2014), and ∼ 0.1 − 0.5 for the multiple impact model, Rufu et al. 2017), while more energetic models, such as the half-Earths model and synestia model, produce nearly pure vapor disks (VMF∼ 0.8 − 1, Nakajima & Stevenson 2014).The VMF of the Moon-forming disk can significantly impact the Moon accretion process; if the initial Moon-forming disk is vapor-poor (liquid-rich), moonlets can quickly form by gravitational instability from a liquid layer in the disk midplane outside the Roche radius (Thompson & Stevenson 1988;Salmon & Canup 2012).These moonlets eventually accrete to the Moon in 10s to 100s of years (Thompson & Stevenson 1988;Salmon & Canup 2012;Lock et al. 2018).In contrast, an initially vaporrich disk needs to cool until liquid droplets emerge before moonlet accretion begins.Growing moonlets in the vapor-rich disk experience strong gas drag from the vapor (Nakajima et al. 2022).The gas drag effect is strongest when the gas and moonlet are coupled to an adequate degree and is sensitive to the moonlet radius.It is strongest when a moonlet radius R p of a few km (Nakajima et al. 2022).In contrast, much smaller moonlets are completely coupled with gas and much larger moonlets are completely decoupled from the gas, experiencing weaker gas drag effect.As a result, ∼km-sized moonlets lose their angular momentum and inspiral onto the Earth within one day (Nakajima et al. 2022), a much shorter timescale than that required for lunar formation.
This same problem was a major challenge to planet formation in the protoplanetary disk (Adachi et al. 1976;Weidenschilling 1977).Here, we quickly review the gas drag problem in the protoplanetary disk.The radial velocity of a particle in a disk is (see Equation (A.6) for the derivation) (Armitage 2010;Takeuchi & Lin 2002), where τ f = Ωt f is the dimensionless stopping time (Armitage 2010), Ω is the Keplerian angular velocity, t f is the friction time (see further descriptions below), v r,g is the gas velocity, η is the pressure gradient parameter described as where c s is the sound speed, and v K is the Keplerian velocity, p g is the gas pressure, and r is the radial distance.v r is largest when τ f = 1.The friction time is the time until the particle and gas reaches the same velocity and is defined as t f = mv rel /F D , where m is the particle mass, v rel is the relative azimuthal velocity between the gas and particle, and F D is the drag force.In the protoplanetary disk, the gas drag for a particle is generally written as , where C D is the gas drag coefficient and ρ g,0 is the initial gas density at the midplane, and R p is the particle radius.The gas drag coefficients are roughly (Armitage 2010) Re < 1 (Stokes regime) 24Re −0.6 , 1 < Re < 800 (Transition regime) 0.44, Re > 800 (Newton regime) (3) where Re is the Reynolds number.Assuming that v rel ∼ ηv K (see Equation A.2), ν ∼ c s λ, where ν is the kinematic viscosity, λ is the mean free path (λ = kBT √ 2πd 2 pg,0 , where k B is the Boltzmann constant, T is the temperature, d is the molecular diameter, p g,0 is the gas pressure at the mid-plane), Re= v rel Rp ν ∼ 4.26, which is in the transition regime, in the protoplanetary disk at 1 AU, assuming T = 280 K, η = 0.002, d = 289 pm, R p = 0.52 m (see discussion below), p g,0 = ρg,0RT mm , and m m = 0.002 kg/mol, where m m is the mean molecular weight and R is the gas constant.Here we also assume the following vertical distribution of the gas density ρ g , where z is the vertical coordinate and H is the gas scale height.This leads to ρ g,0 = Σg √ 2πH .We assume Σ g = 17000 kgm −2 and use the relationships of 2 .The orbital angular frequency is Ω = GM/r 3 , where M is the stellar mass, and r is the distance from the star (for the Moon-forming disk, M is the Earth mass and r is the distance from Earth).These parameters are taken from previous work (Carrera et al. 2015) (other values of β have been discussed, such as β = 3 7 , Chiang & Youdin 2010).
The dimensionless stopping time becomes τ f = 1 when R p = 0.52 m.The particle density ρ p = 3000 kgm −3 is assumed.The residence time of the particle at 1 AU is 1 AU/v r = 75 years, where v r is the radial fall velocity of the particle (see A.7.For simplicity, the radial gas velocity v r,g = 0 is assumed).Thus, an approximately meter-sized particle falls toward the Sun at 1 AU within 80 years, a timescale much shorter than the planet formation timescale (several-tens of Myr).This is the so-called "meter barrier" problem, which was a major issue to explain planetary growth in a classical planet formation model (Adachi et al. 1976;Weidenschilling 1977).
In contrast, for the Moon-forming disk, τ f = 1 (Re > 10 10 ) is achieved when R p = 1.8 km, assuming ρ g,0 = 40 kg/m 3 , ρ p = 3000 kg/m 3 , r = 3R ⊕ , η = 0.04, T = 5200 K (see Section 2.3 for justifications for these parameters), and d = 300 pm.The residence time, r/v r , is approximately 1 day (1.21 day), which is much shorter than the Moon-formation timescale of several 10s of years to 100 years.This formation timescale is ultimately determined by the radiative cooling timescale, but it is model dependent.Here we provide a very simple estimate; the time scale for radiative cooling is 4πr 2 σT 4 ph ∼ 10 years (this is also consistent with numerical work, Lock et al. 2018), where M disk is the disk mass, L is the latent heat, C p is the specific heat, ∆T is the temperature change over time, σ is the Stefan-Boltzmann constant, and T ph is the photosphere temperature.Here, M disk ∼ 0.015M ⊕ , L = 1.2 × 10 7 J/kg (Melosh 2007), C p = 10 3 J/K/kg, ∆T = 2000 K, T ph = 2000 K (Thompson & Stevenson 1988) are assumed.However, the actual Moon-formation timescale can be longer than this for several reasons, including additional heating due to viscous spreading (Thompson & Stevenson 1988;Charnoz & Michaut 2015), and radial material transport efficiency (e.g.Salmon & Canup 2012).Short Moon-formation timescale has been proposed (e.g.Mullen & Gammie 2020;Kegerreis et al. 2022), but the typical Moon-formation timescale has been estimated to range in 10 − 10 2 years.The vertical settling velocity of condensing particles ρpRpΩ 2 z ρg (Armitage 2010).Assuming z ∼ H, the settling time for a particle with a radius of 2 km is 0.22 days.Determining the collision history of moonlets require conducting orbital dynamics simulations, but here we provide a rough estimate.A rough estimate of the mass doubling time of the largest moonlet in the Moon-forming disk is typically ∼ 1 day (Salmon & Canup 2012).If a 2 km-sized moonlet mass doubles in a day, the radius change is very small (2km×2 1/3 = 2.5 km) and the gas drag effect remains strong on such a moonlet and this change does not prevent the moonlet from falling into Earth.However, the actual collision time can vary depending on the local concentration of particles, which needs a detailed future study.
This gas drag effect is a problem with forming the Moon from an initially vapor-rich disk.The Moon can still form after most of the vapor condenses, but by that time a significant portion of the disk mass could be lost (Ida et al. 2020;Nakajima et al. 2022), which would fail to form a large moon from an initially vapor-rich disk (here we use a "large" moon when its mass is ∼ 1 wt% or larger of the host planet).If this is the case, an initially vapor-rich disk may not be capable of forming a large Moon (Nakajima et al. 2022).In contrast, the gas drag effect is weak for vapor-poor disks, which are generated by less energetic models, such as the canonical and multiple impact models.

Streaming instability in the general impact-induced moon-forming disk
A potential solution to this vapor drag problem is forming a large moonlet very quickly (much larger than 2 km), so that the moonlet would not experience strong gas drag.This is the accepted solution for the gas drag problem in the protoplanetary disk.The proposed mechanism is the streaming instability (Youdin & Goodman 2005;Johansen et al. 2007), where particles spontaneously concentrate in the disk, gravitationally collapsing and forming a large clump (∼ 100 km in size, Johansen et al. 2015).If this mechanism works for the Moon-forming disk, an initially vapor-rich disk may be able to form the Moon despite the gas drag issue.If this mechanism turns out not to work for the Moon-forming disk, it is an interesting finding as well, given that a Moon-forming disk is often treated as a miniature analogue of the protoplanetary disk.Understanding what makes these two disks differ would deepen our understanding of planet and satellite formation processes.
Moreover, knowing whether the streaming instability can affect moon formation processes informs our understanding of moon formation in the solar and extra-solar systems.Moon formation in an impact-induced disk is common in the solar system (e.g.Martian moons, Craddock 2011, Uranian moons, Slattery et al. 1992, andPluto-Charon Canup 2005).While there are no confirmed exomoons (moons around exoplanets) to date (e.g.Cassese & Kipping 2022), impact-induced exomoons should be common because impacts are a common process during planet formation (Nakajima et al. 2022) and because impacts in extrasolar system may have been already observed (e.g.Meng et al. 2014;Bonomo et al. 2019;Thompson et al. 2019;Kenworthy et al. 2023).If streaming instability operates in these disks, it can affect what types of planets can host exomoons, which can be compared with future exomoon observations.
It should be also noted that other instability, such as secular gravitational instability (GI) (e.g.Youdin 2011;Takahashi & Inutsuka 2014;Tominaga et al. 2018) and two-component viscous GI (TVGI) (Tominaga et al. 2019), have been discussed as a mechanism for planetesimal and dust ring formation.The GI occurs when dust-gas interaction reduces the rotational support of the rotating disk, which leads to dust concentration.The TVGI is an instability cased by dust-gas friction and turbulent gas viscosity.These instabilities can lead to clump formation even if the disk is self-gravitationally stable.Implications of these instabilities are beyond the scope of this paper.

Motivation of this work
The goal of this work is to investigate whether the streaming instability can form moonlets that are large enough to avoid strong gas drag from the vapor-rich disk.The result would constrain the Moon-formation model as well as general impact-induced models in the solar and extrasolar systems.We conduct hydrodynamic simulations using the code Athena (Stone et al. 2008;Bai & Stone 2010a).We first conduct 2D simulations to identify the section of parameter space that leads to streaming instability (SI).Subsequently, we conduct 3D simulations with self-gravity to identify the size distribution of moonlets.Lastly, we identify the lifetime of the moonlets formed by streaming instability to investigate whether it is possible to form a large moon from an initially vapor-rich disk.For the general impact-induced disks, we consider "rocky" and "icy" disks, where these disks form by collisions between rocky planets (with silicate mantles and iron cores) and between icy planets (with water ice mantles and silicate cores), respectively.

Athena
We use the Athena hydrodynamics code, which solves the equations of hydrodynamics using a second-order accurate Godunov flux-conservative approach (Stone et al. 2008).We use the configuration of Athena that couples the dimensionally-unsplit corner transport upwind method (Colella 1990) to the third-order in space piecewise parabolic method by Colella & Woodward (Colella & Woodward 1984) and calculates the numerical fluxes using the HLLC Riemann solver (Toro 1999).We also integrate the equations of motion for particles following Bai & Stone (Bai & Stone 2010a) and include particle self-gravity for 3D simulations following the particlemesh approach described in previous work (Simon et al. 2016).Orbital advection is taken into account following previous work (Bai & Stone 2010a,b).Our setup is the local shearing box approximation in which a small patch of the disk is corotating with the disk at the Keplerian velocity (Stone & Gardiner 2010).The local Cartesian frame is defined as (x, z) for 2D and (x, y, z) for 3D simulations, with x as the radial coordinate from the planet and z is parallel to the planetary spin axis, and y is in the direction of orbital rotation.
In the Athena code, we solve the following equations for our simulations (Bai & Stone 2010a;Simon et al. 2016;Li et al. 2019): (8) The first two equations specify mass and momentum conservation for the gas, respectively, while the third equation represents the motion of a particle i coupled with the gas.Here, u is the velocity of the gas and I is the identity matrix.v is the mass-weighted averaged particle velocity in the fluid element, assuming that particles can be treated as fluid (Bai & Stone 2010a).The terms of the right hand side of Equation ( 7) are radial tidal forces (gravity and centrifugal force), vertical gravity, and the Coriolis force, and the feedback from the particle to the gas.In Equation ( 8), the first term on the right hand side describes a constant radial force due to gas drag.v i is the particle velocity, x and ẑ represent the unit vectors in the x and z directions, x i and z i are the values of x and z for the particle i. u is the gas velocity interpolated from the grid cell centers to the location of the particle.The second, third, and forth terms are radial and vertical tidal forces and the Coriolis force.a g is the acceleration due to self-gravity, which is considered only in 3D simulations.a g = −∇Φ p , where Φ p is the gravitational potential and satisfies Poisson's equation, ∇ 2 Φ p = 4πGρ p .To reduce computational time, the particles are organized into "superparticles", each representing a cluster of individual particles of the same size.In the code units, we normalize Ω = c s = H = ρ g,0 = 1.The gas and particle initial distributions are described as Equation ( 4) and respectively, where H p is the scale height of particles and is set to 0.02H (Bai & Stone 2010a) and Σ p is the particle surface density.The system uses an isothermal equation of state P = ρ g c 2 s , and the particles are distributed uniformly in the x and y direction and normally in the z direction (Equation 9).The computational domains are 0.2H × 0.2H in 2D simulations and 0.2H × 0.2H × 0.2H in 3D simulations, with all-periodic boundary conditions.The resolution for 2D is 512 × 512 and each grid cell has 1 particle (512 × 512 = 262, 144 particles) .The resolution for 3D is 10 × 128 3 (= 21 million particles in total).Previous work shows that the resolution of 128 3 can produce the large clump size distribution well compared to those with 256 3 and 512 3 (Simon et al. 2016) and therefore this 3D resolution is sufficient to resolve large clumps, which are the focus of this work.The initial particle size is assumed to be constant.Variable initial particle sizes could affect the growth speed (Krapp et al. 2019) and the concentrations of particles depending on their sizes (Yang & Zhu 2021), but it is not known to affect the largest clump size formed by streaming instability.

Clump detection in 2D and 3D
After we conduct 2D simulations, we identify filaments forming in the disk in order to constrain the parameter space favorable for the streaming instability.We use the Kolmogorov-Smirnov (KS) method, which is based on previous work (Carrera et al. 2015).For each time step in a given time window (25/Ω in our case), we compute the particle surface density Σ p (x) = 0.1 −0.1 ρ p (x, z) dz and average it over the time window to give ⟨Σ p ⟩ t .We then sort the values in ⟨Σ p ⟩ t from highest to lowest, and compute the cumulative distribution of this sorted dataset.Finally, we use the Kolmogorov-Smirnov test to output a p-value that measures the likelihood that the underlying cumulative distribution was linear: where D is the maximum distance between the data and linear cumulative distributions, and n is the number of data points.The higher the p is, the more homogeneous the system is and therefore filament formation is unlikely, whereas small p indicates filament formation is more likely.Here we assume that filament formation is very likely when p < 0.10, the same as in previous work (Carrera et al. 2015).
For self-gravitating clump detection in our 3D simulations, we use PLanetesimal ANalyzer (PLAN) (Li et al. 2019).This is a tool to identify self gravitating clumps specifically made for Athena output.The density of each particle is assessed based on nearest particles.Particles with densities higher than a threshold are associated with neighbouring dense particles until a density peak is achieved.In contrast, particle groups with a saddle point less than a threshold remain separated.Detailed descriptions are found in previous work (Li et al. 2019).

Model Parameters
The main parameters for the Athena simulations are (1) the dimensionless stopping time τ f (the Newton regime, see Section 1.2 and Nakajima et al. 2022), (2) the normalized pressure parameter ∆ = ηv K /c s , (3) the ratio of the particle surface density to the gas surface density Z = Σ p /Σ g (VMF= 1 1+Z ), and (4) the normalized gravity G ≡ 4πGρ g,0 /Ω 2 for 3D simulations, where G is the gravitational constant.A small value of τ f indicates a small particle radius R p , where the particle is well coupled with the gas, whereas a large value of τ f corresponds to a large value of R p , which is more decoupled from the gas.A large value of ∆ corresponds to a large pressure gradient and quicker radial infall.G represents the strength of self-gravity.The parameter space we are exploring for our 2D simulations is τ f = 10 −3 , 10 −2 , 10 −1 , 10 0 , 10 1 , ∆ = 0.1, 0.2, 0.3, 0.4, 0.5 and Z = 0.05, 0.1.For 3D simulations, we use Z = 0.1, 0.3 and G = 0.1788 and 0.5898, where the former G value corresponds to slightly cooler temperature (4700 K) while the latter corresponds to hotter temperature (5200 K).Here, we justify the choice of these parameters.This range of τ f corresponds to the particle radius of 2 m and 20 km (see Equation 5).The global disk structures have been calculated based on hydrodynamic simulations in previous work.The overall disk mass of the Moon-forming disk is typically a few percent of Earth, depending on the impact model (M D /M L =1.35-2.80,where M D is the disk mass and M L is the lunar mass, Nakajima & Stevenson 2014).The mid-plane disk temperature ranges from 3000 K to 7000 K and the radial range of the disk is ∼1-8 R ⊕ (see Figure 5 in Nakajima & Stevenson 2014).The disk tem-perature is ∼ 4000−5500 K at r = 3R ⊕ .For the general rocky and icy impact-induced disks, the disk temperature can vary, but typically in the range of thousands of K for vapor-rich disks (Nakajima et al. 2022).The pressure gradient can vary, but the typical value of η is ∼ 0.02 − 0.06 based on impact simulations (Nakajima et al. 2022).In the Moon-forming disk, the value of Z increases as the disk cools (on the timescale of 10 − 10 2 years).In other words, Z can be zero initially in an energetic Moon-forming impact model (e.g.Lock et al. 2018) and eventually becomes infinity as the disk materials condense and the gas disappears.Thus, picking a value of Z means that we are seeing physics at a specific time.Since we are primarily interested in high-vapor disks (i.e. an early phase of the disk), we explore the mass ratio range Z ∈ [0.05, 0.1] for our 2D simulations and Z ∈ [0.1, 0.3] for our 3D simulations.We use the large value of Z = 0.3 as a sensitivity test.Higher values of Z can be achieved as the disk cools, but we focus on small values of Z(≤ 0.3) for several reasons.First, at a larger value of Z(> 0.3), the conventional gravitational instability in the liquid part of the disk can occur to form moonlets.At Z = 0.3, VMF= 1 Z+1 = 0.76 and the total (Σ p +Σ g ) surface density at r ∼ 3R ⊕ is ∼ 10 8 kg m −2 in energetic models, which means that Σ p = 0.76 × 10 8 kg m −2 = 7.6 × 10 7 kg m −2 .Previous work suggests that when the liquid (melt) layer's thickness reaches 5-10 km (or equivalent of a few 10 7 kg m −2 ), gravitational instability can happen in the melt layer (Machida & Abe 2004).Whether streaming instability occurs at the same time, or whether streaming instability affects the gravitational instability in the Moon-forming disk are unknown and have not been explored in previous studies.Under these circumstances, the important of streaming instability becomes unclear.Additionally, these streaming instability simulations have been conducted at low Z values (≤ 0.1) in previous work to reproduce conditions in the protoplanetary disks (e.g.Abod et al. 2019;Li & Youdin 2021) and simulations with high Z(> 0.3) values have not been fully tested.For these reasons, we focus on relatively small values of Z in this study.
For the set of Athena simulations, we focus on reproducing two sets of the Moon-forming disk thermal profile; assuming M = M ⊕ and r = 3R ⊕ , where M ⊕ is the Earth mass and R ⊕ is the Earth radius, the gas pressure at the midplane around 3R ⊕ is ≈ 12 MPa, T = 4700 K, and m m = 30 g/mol.This makes the density ρ g,0 = 12.13 kg m −3 for an ideal gas, c s = 1140 m/s, v K = 4562 m/s, Ω = 2.38 × 10 −4 s −1 , and the scale height H = c s /Ω = 4.78 × 10 6 m, and G = 0.1788.For the higher-temperature scenario, T = 5200 K, ρ g,0 = 40.01kg m −3 , c s = 1200 m/s, H = 5.03×10 6 , and G = 0.5898.These temperatures are motivated by previous hydrodynamic calculations of the Moon-forming disk formation (Nakajima & Stevenson 2014) and other parameters are calculated based on an equation of state of dunite assuming that the vapor and liquid phases are in equilibrium (MANEOS, Thompson & Lauson 1972;Melosh 2007).Since we are assuming the disk is in the liquid-vapor equilibrium, higher disk temperature leads to higher vapor density at the mid-plane.For example, the gas density at the liquid-vapor equilibrium is 150 kg/m 3 at 6000 K and 1.8 kg/m 3 at 4000K for dunite, according to MANEOS (Thompson & Lauson 1972;Melosh 2007).Assuming η ∼ 0.02 − 0.06 for the Moon-forming disk, ∆ ∼ 0.1 − 0.2.However, higher values are possible (Nakajima et al. 2022), and therefore we explore the range of ∆ ∼ 0.1 − 0.5.
The parameter values of ∆ and G are significantly different from values in the protoplanetary disk, where the typical values used are ∆ ∼ 10 −3 −10 −2 and G ∼ 0.05 in the protoplanetary disk (Carrera et al. 2015;Simon et al. 2016).We use a fixed value of Z in hydrodynamic calculations because the condensation timescale (∼years) is longer than the simulation timescale (we run 2D simulations for ∼ 100 orbits, which correspond to 123 hours.A steady state is reached by this time).We also assume that the disk does not evolve on this short timescale.

2D Athena simulations
Figure 1A shows four examples of space time plots of our 2D simulations, where (1) (τ f , ∆) = 10 −3 , 0.3, (2) 10 −2 , 0.1, (3) 10 −1 , 0.2, and (4) 10 −1 , 0.2.Z = 0.1 for the four cases.These simulations represent four characteristic regimes.The horizontal axis is x normalized by the scale height H and the vertical axis is the number of orbits.The color shows particle concentration.In case (1), τ f is small, which means that gas and particles are well coupled and this does not lead to filament formation.In case (2), after ∼ 20 orbits, filaments start forming and these filaments are stable during the rest of the simulations, which indicate streaming instability occurs with this chosen parameter.A radially shearing periodic boundary condition is used in the x direction, and therefore the filament that appears at x/H = −0.1 at ∼ 100 orbits reappears at x/H = 0.1.Two distinct filaments form in this simulation.In case (3), filaments form, but they are not stable because their radial movements are high due to the high ∆ value.Thus, this is not an ideal parameter space for SI.Similar behaviors have been seen in simulations for the protoplanetary disk with other parameter combinations (e.g.∆ = 0.05, τ f = 3, Z = 0.005, Carrera et al. 2015).In case (4), filaments are not as clearly defined as those in case ( 2), but a coherent filament structure is found after several orbits.Thus, this is also considered as an SI regime.
Figure 1B shows the summary of our 2D simulations for Z = 0.1 (top) and for Z = 0.05 (bottom).We use the Kolmogorov-Smirnov test (KS) to identify the extent of particle concentration (see Section 2.2).Concentration is measured by the p value and when the value of p is small (p < 0.1), streaming instability is likely (Carrera et al. 2015).The colorbar indicates the p value.At Z = 0.1, p < 0.1 is achieved when ∆ = 0.1 and 10 −3 ≤ τ f ≤ 10 0 .At large ∆ values, radial motions are large and stable filaments do not form, as discussed above.At Z = 0.05, this trend remains the same, but the streaming instability regime is smaller (10 −2 ≤ τ f ≤ 10 0 ) than that of Z = 0.1.This is because larger Z leads to more particles in the disk and therefore to more filament formation (Carrera et al. 2015).Thus, our 2D simulations show that the most filaments form at relatively small ∆ (∆ = 0.1) and with both Z = 0.1 and Z = 0.05, but the higher Z has more favorable conditions.This general trend is consistent with previous work in the protoplanetary disk (e.g.Carrera et al. 2015), which investigate smaller ∆ values for the protoplanetary disk (∆ ≤ 0.05).

3D Athena simulations
Now that we identify the parameter space for the streaming instability (∆ = 0.1, Z = 0.1), we perform 3D simulations with self-gravity to estimate the size distributions of SI-induced moonlets.Figure 2 shows snapshots of one of our 3D simulations (τ f = 1, G = 0.5898, Z = 0.1 and ∆ = 0.1) at four different times (t = 0.32, 2.87, 3.39, 3.55), where t is the number of orbits.The self-gravity is not initially included and is turned on at t = 3.18.This is a general practice to avoid artificial clumping before streaming instability takes place (Simon et al. 2016).The color shows the particle density.The circles indicate the Hill spheres of self-bound clumps, which are identified using PLAN (Li et al. 2019).At t = 0.32, no clump is identified, but by t = 2.87, the streaming instability fully develops and reaches its steady state.After self gravity is turned on at t = 3.18, self-gravitating clumps form right away (t = 3.39).This general behavior is similar to previous work on the streaming instability in the protoplanetary disk (Abod et al. 2019), but the time it takes to form clumps by streaming instability is shorter, probably because of the large Z value (Simon et al. 2022) compared to those of the protoplanetary disk.
Figure 3A shows the cumulative distribution of moonlet mass M p , normalized by the characteristic self- gravitating mass M G , defined as (Abod et al. 2019), where λ G is an instability wavelength, which originates from the Toomre dispersion relation, equating the tidal and gravitational forces (Abod et al. 2019).The parameters are τ f = 0.1, 1, G = 0.1788, 0.5898 and Z = 0.1, 0.3.For the G = 0.1787 cases, M G = 5.19 × 10 18 kg and for the G = 0.5898 cases, M G = 2.17 × 10 20 kg.For all the 3D simulations, ∆ = 0.1 are assumed.The solid lines indicate the moonlet mass distribution shortly after the onset of the self gravity (t = 3.66 for G = 0.1788, τ f = 0.1 and t = 3.34 for all the other cases).The dashed lines indicate the same parameter after an additional time (t + ∆t, where ∆t = 1/Ω = 0.16 orbit except the Z = 0.3 case where ∆t = 1 2Ω .The reason of the shorter ∆t for the higher Z case is that the streaming instability develops faster for higher Z values, Simon et al. 2022).The maximum M p /M G = 0.254 ( G = 0.1788, τ f = 1 at t = 3.50).This is broadly consistent with previous studies that suggest that the maximum clump masses formed by the streaming instability are characterized by M G , ranging 10 −1 − 10 1 M G (Abod et al. 2019;Li et al. 2019).Our result lies on the lower end of this mass range.This indicates that the moonlet mass distribution based on our numerical simulations is consistent with the analytical mass model generated for the protoplanetary disk.Thus, this also indicates that the processes of streaming instability are similar between the Moon-forming disk and protoplanetary disk despite the different input parameters.We also find that higher Z does not necessarily lead to higher M p .
Here we consider the best case scenario of forming a large moonlet.Assuming that the density of the clump is 3000 kg m −3 , then M G = 2.17 × 10 20 kg is equivalent to 258 km in radius and 0.254 M G = 163 km.Based on equation (1), the residence time of 258 km and 163 km moonlets are 92 and 58 days, respectively, assuming v r,g = 0.These timescales are much shorter than the Moon-formation timescale (10s-100 years, Thomp- At first (t = 0.32), no concentration of particles is observed, but streaming instability clearly develops by t = 2.87.After self-gravity is turned on at t = 3.18, moonlets form by gravitational instability (t = 3.39, 3.55).The self-gravitating clumps are identified using PLAN (see the main text).
son & Stevenson 1988).Thus, even though the streaming instability can occur in an initially vapor-rich Moonforming disks, it does not help increasing the residence time of moonlets.

Exomoon formation by streaming instability
The streaming instability likely occurs in impactinduced moon-forming disks in extrasolar systems.The pressure gradient parameter η (∼ 0.02 − 0.06) is similar regardless of the composition of the disk (Nakajima et al. 2022).Figure 3B shows the disk residence time for a clump with mass M G in a moon-forming disk as a function of the planetary mass (M planet = 1 − 6M ⊕ )."Rocky Planet" corresponds to disks formed by a collision between terrestrial planets while "Icy Planet" indicates disks formed by a collision between icy planets whose mantles are made of water ice (70 wt%) and cores are made of iron (30 wt%).This is produced by calcu-lating r/v r (see Equation 1), where for rocky planets we use ρ g,0 = 40.01kg m −3 , T = 5200 K, ρ p = 3000 kg m −3 , r = 3R planet where R planet is the planetary radius, and for icy planets we use ρ g,0 = 10.0 kg m −3 , T = 2000 K, ρ p = 1000 kg m −3 , r = 3R planet .These are the temperatures when the impact-induced disks reach complete vaporization based on impact simulations (Nakajima et al. 2022).Higher temperature is possible, but the disk would need to cool down to this temperature to form particles (dust), and therefore this is the most relevant temperature to assess the effect of streaming instability.The thermal state of the disk formed by icy planetary collisions are estimated based on an equation of state of water (Senft & Stewart 2008).Additionally, R planet = R ⊕ (M planet /M ⊕ ) 1/3.7 for rocky planets and R planet = 1.2R ⊕ (M planet /M ⊕ ) 1/3 for icy planets are assumed (Kipping et al. 2013;Mordasini et    2012).At M planet = 1 − 6M ⊕ , these parameters make G = 0.42 − 0.67 for rocky planets and G = 0.25 for icy planets, which are similar to the ranges covered in our hydrodynamic simulations (Section 2.3).
In both icy and rocky planet cases, the disk residence time is a few 10s of days to several months, which are short compared to the satellite formation timescale (10s-100s years, Nakajima et al. 2022).Another effect is that as the disk cools, ρ g,0 and H decrease, which means that M G decreases over time in both the icy and rocky disks.This means that SI-induced clumps tend to become smaller as time progresses.These effects indicate that the streaming instability likely plays a limited role in impact-induced moon-forming disks.Figure 4 shows a schematic view of our hypothesis.An energetic impact would generate a vapor-rich disk (the disk is made of silicate vapor for the Moon or rocky planets and water vapor for icy planets).Over time, the disk cools by radiation and small droplets (< cm) emerge.These small droplets would grow by accretion and by streaming instability.However, once these moonlets reach 100 m -100 km in size, gas drag from the vapor is so strong that they fall onto the planet on a short time scale (days-weeks).This continues to occur until the vapor mass fraction of the disk decreases so that the gas drag effect is no longer strong.Once this condition is reached, a liquid layer emerges and moonlets can stay in the disk.However, by this time a significant disk mass could have been lost (> 80 wt %) (Ida et al. 2020;Nakajima et al. 2022).For this reason, only a small moon (or moons) could form from an initially vaporrich disk.Thus, we suggest that an initially vapor-rich disk is not suitable for forming our Moon, and our result supports the hypothesis that the Moon formed from an initially vapor-poor disk, including the canonical model where the proto-Earth was hit by a Mars-sized impactor (Canup & Asphaug 2001).
The isotopic composition of the Moon may constrain whether streaming instability played a role during the Moon formation.Some may argue that the observational fact that the Moon is depleted in volatiles may be explained by accreting moonlets formed by streaming instability before volatiles accreted onto the Moon.However, this process would make the Moon isotopically light if these isotopes experience kinetic fractionation (Dauphas et al. 2015(Dauphas et al. , 2022)), which is inconsistent with the observation of enrichment of heavy isotopes in the Moon (such as K, Wang & Jacobsen 2016).The lunar isotopes would be heavier if they experience equilibrium fractionation, but the equilibrium fractionation would not produce the observed isotopic fractionations at the high disk temperature condition.Alternatively, some volatiles could be lost from the lunar magma ocean under an equilibrium condition (Charnoz et al. 2021), but the efficiency of volatile loss is not fully known (Dauphas et al. 2022).

Streaming instability in general impact-induced disks
Our model supports the previous work that suggests that relatively large rocky (> 6M ⊕ , > 1.6R ⊕ ) and icy (> 1M ⊕ , > 1.3R ⊕ ) planets cannot form impactinduced moons that are large compared to the host planets (Nakajima et al. 2022).Larger planets than those thresholds generate completely vapor disks, because the kinetic energy involved in an impact scales with the planetary mass.Thus, these large planets are not capable of forming moons that are large compared to their host planets.Moons can form by mechanisms other than impacts, such as formation in a circumplanetary disk and gravitational capture, but these moons tend to be small compared to the sizes of their host planets (the predicted moon to planet mass ratio is ∼ 10 −4 for moons formed by circumplanetary disk Canup & Ward 2006 and gravitationally captured moons are small in the solar system, Agnor & Hamilton 2006).Thus, fractionally large moons compared to the host planet sizes, which are observationally favorable, likely form by impact.So far, no exomoon has been confirmed despite extensive searches, but future observations, especially with James Webb Space Telescope (Cassese & Kipping 2022), may be able to find exomoons and test this theoretical hypothesis.

Comparison between the Moon-forming disk and protoplanetary disk
It is certainly intriguing that the streaming instability is able to solve the gas drag problem in the protoplanetary disk, but not in the Moon-forming disk.In both scenarios, the clump sizes formed by the streaming instability happen to be similar (M G ∼100 km, see Johansen et al. 2015, for the protoplanetary disk and ∼100 km for the protolunar disk, respectively).This size is ∼ 10 5 times larger than the size of the particle (0.52-1 m) (Adachi et al. 1976) that experiences the strongest gas drag (τ f = 1 and v r ∼ −ηv K ).Therefore, once particles become ∼ 100 km-sized planetesimals by streaming instability, the radial velocity of the planetesimals decrease drastically (v r = − 2ηvK τ f for large τ f , see Equation A.7), which helps planetesimal growth.In contrast, the largest moonlet size (100 km) formed by streaming instability in the Moon-forming disk is only 50 times larger than the particle size that experiences the strongest gas drag (2km).This does not result in a large τ f change or v r change, and therefore streaming instability does not effectively help moon formation.Therefore, streaming instability is an effective mechanism to skip the "1 m barrier" in the protoplanetary disk, whereas it is not an effective mechanism to skip the "km barrier" in the Moon-forming disk.

Streaming instability in the circumplanetary disk around a gas giant
The streaming instability may occur during moon formation in various satellite systems around gas giants.The Galilean moons around Jupiter and Titan around Saturn likely formed from their circumplanetary disks.In these disks, the particle-to-gas ratio Z is typically 10 −4 − 10 −2 , which was thought to be too small for streaming instability to occur (Carrera et al. 2015;Yang et al. 2017;Shibaike et al. 2017).However, recent SI calculations show that SI can occur at as low as Z = 4 × 10 −3 (at ∆ = 0.05 and τ f = 0.3, Li & Youdin 2021, and the required Z depends on the disk conditions, e.g.Carrera et al. 2015;Yang et al. 2017;Sekiya & Onishi 2018) which may indicate that SI can be potentially important in circumplanetary disks around gas giants if the dust-gas ratio of the disk is relatively high.The velocity difference between the gas and particles v rel is (Canup & Ward 2002) Assuming cs rΩ = cs v k ∼ 0.1 (Canup & Ward 2002), this yields η ∼ 0.01 and ∆ = ηv k cs ∼ 0.1.If Z is sufficiently large (at least Z > 4 × 10 −3 Li & Youdin 2021), it is possible that the streaming instability takes place in a circumplanetary disk around a gas giant.A rough estimate for the size of a moonlet in a circumplanetary disk is M G = 1.4×10 16 kg (Equation 12), which is ∼ 10.3 km in radius if it is a rocky moonlet.Here, we are assuming Z = 10 −2 , Σ p = 3 × 10 4 kg m −2 , c s = 1000 m s −1 , H ∼ c s /Ω, ρ g = Σ g /(H √ 2π) (Canup & Ward 2002).These moonlets are relatively small compared to large moons around gas giants (e.g. the Galiean moon masses range 10 22 − 10 23 kg) and it is unclear if they would significantly impact these moon-formation process.Further research is needed to understand its effect on the satellite formation in a circumplanetary disk around gas giants.

Model limitations
There are several model limitations that need to be addressed in our future work.First, the effect of the Roche radius (a R ∼ 3R ⊕ ) needs to be taken into account to understand the Moon formation.An inspiraling moonlet would not directly reach the Earth, but would be tidally disrupted near the Roche radius, where it then might be incorporated into an accretion disk (e.g., Salmon & Canup 2014).Secondly, our model presented here does not take into account the evolution of the disk in detail, which is important to identify the final mass and composition of the resulting moon.As the disk spreads out, it is likely that ∆ decreases, which would slow down the radial drift of moonlets.This gas drag effect would disappear once the local vapor condenses, which would occur at the outer part of the disk first, since that part of the disk can cool efficiently due to its large surface area.Given that such moonlets that directly form from the disk would not have time to lose volatiles in the disk phase and therefore this model requires that the Moon lost its volatiles before or after the disk phase, for example during the lunar magma ocean phase (Charnoz et al. 2021) (see further discussion in Section 4.1).Nevertheless, such scenario may be possible in moon-forming disks in the solar and extrasolar systems.The effects of the Roche radius and disk evolution will be addressed in our future work.

CONCLUSION
In conclusion, we show that the streaming instability can form self-gravitating clumps (∼ 10 2 km) in a vaporrich moon-forming disk generated by a giant impact, but the sizes of the clumps it generates are not large enough to avoid inspiraling due to the strong gas drag.This is a major difference from the protoplanetary disk, where the streaming instability can efficiently form large clumps to avoid the strong gas drag effect.As a result, growing moonlets in an initially vapor-rich moon-forming disk continue to fall onto the planet once they reach the sizes of 100 m -100 km.These moonlets could grow further once the disk cools enough and the vapor mass fraction of the disk becomes small.However, by this time a significant amount of the disk mass is lost, and the remaining disk could make only a small moon.This result is applicable to impact-induced moon-forming disks in the solar system and beyond; we find that the streaming instability is not an efficient mechanism to form a large moon from an impact-induced vapor-rich disks in general.As a result, we support previous work that suggests that fractionally-large moons compared to their host planets form from vapor-poor disks; the ideal planetary radii that host fractionally large moons are 1.3 − 1.6R ⊕ (Nakajima et al. 2022) given that rocky or icy planets larger than these sizes would likely produce completely vapor disks, which are not capable of forming large moons (Nakajima et al. 2022).The streaming instability may take place in circumplanetary disks, but their effect on the moon-formation process needs further investigation.

Figure 2 .
Figure 2. Snapshots of a 3D simulation (τ f = 1, ∆ = 0.1, Z = 0.1, G = 0.5898).The horizontal and vertical axes are x and y, normalized by the scale height H.The color indicates the particle density, normalized by the average density.t indicates the number of orbits.The circles represent the Hill radius of each moonlet formed by streaming instability.At first (t = 0.32), no concentration of particles is observed, but streaming instability clearly develops by t = 2.87.After self-gravity is turned on at t = 3.18, moonlets form by gravitational instability (t = 3.39, 3.55).The self-gravitating clumps are identified using PLAN (see the main text). al.

Figure 3 .
Figure 3. (A) Cumulative mass distribution of moonlets formed by streaming instability at ∆ = 0.1, Z = 0.1, 0.3.The purple, dark blue, light blue, dark yellow, and light yellow lines indicate parameter values of ( G, τ f , Z) = (0.5898, 0.1, 0.1), (0.1788, 0.1, 0.1), (0.5898, 1, 0.1), (0.1788, 1, 0.1), (0.1788, 0.1, 0.3).The solid and dashed lines represent the mass distribution at different times (see Section 3.2).The horizontal axis indicates the moonlet mass (Mp) normalized by MG and the vertical axis indicates the number of moonlets whose masses are larger than the given moonlet mass.(B) Residence time of a moonlet whose mass is MG formed by streaming instability as a function of the planetary mass M planet normalized by the Earth mass at the Roche radii.The blue solid and yellow dash-dot lines represent disks formed by collisions between icy planets and rocky planets, respectively.

Figure 4 .
Figure 4. Schematic view of the Moon-formation from an initially vapor-rich disk.
Streaming instability in the Moon-forming disk