How Long-lived Grains Dominate the Shape of the Zodiacal Cloud

Grain–grain collisions shape the three-dimensional size and velocity distribution of the inner Zodiacal Cloud and the impact rates of dust on inner planets, yet they remain the least understood sink of zodiacal dust grains. For the first time, we combine the collisional grooming method combined with a dynamical meteoroid model of Jupiter-family comets (JFCs) that covers 4 orders of magnitude in particle diameter to investigate the consequences of grain–grain collisions in the inner Zodiacal Cloud. We compare this model to a suite of observational constraints from meteor radars, the Infrared Astronomical Satellite, mass fluxes at Earth, and inner solar probes, and use it to derive the population and collisional strength parameters for the JFC dust cloud. We derive a critical specific energy of QD*=5×105±4×105Rmet−0.24 J kg−1 for particles from JFC particles, making them 2–3 orders of magnitude more resistant to collisions than previously assumed. We find that the differential power-law size index −4.2 ± 0.1 for particles generated by JFCs provides a good match to observed data. Our model provides a good match to the mass-production rates derived from the Parker Solar Probe observations and their scaling with the heliocentric distance. The higher resistance to collisions of dust particles might have strong implications to models of collisions in solar and exosolar dust clouds. The migration via Poynting–Robertson drag might be more important for denser clouds, the mass-production rates of astrophysical debris disks might be overestimated, and the mass of the source populations might be underestimated. Our models and code are freely available online.


INTRODUCTION
Despite its tenuous appearance (we can see stars and galaxies from our planet), solar system dust reveals its presence via collisions.Many spacecraft have been directly affected or damaged by particle impacts (Drolshagen & Moorhead 2019;Szalay et al. 2020;Malaspina et al. 2022;Williams et al. 2017) and consequential plasma generation (Close et al. 2010).Surfaces of airless bodies are continuously eroded, a process thought to produce tenuous extraterrestrial exospheres (Killen & Hahn 2015;Pokorný et al. 2019Pokorný et al. , 2021)).
The mutual collisions of dust and meteoroids have been studied for decades, sometimes controversially.Some studies predict that the collisions are a main driver of dust/meteoroid removal in the inner solar system (e.g., Grun et al. 1985;Mann et al. 2004), whereas other studies show that collisions cannot be as frequent as previously thought if models are to match the shape and distribution of the Zodiacal Cloud (e.g., Nesvorný et al. 2011a;Pokorný et al. 2014;Soja et al. 2019).
Before the advent of high-performance computing, mutual meteoroid collisions were examined analytically.Such analytic treatments have succeeded in reproducing several observed meteoroid-related phenomena.For instance, Wiegert et al. (2009) implemented a hybrid of Grun et al. (1985) and Steel & Elford (1986) collisional lifetimes to support a dynamical model that was able to model certain aspects of the orbital distribution of radar meteors at Earth using cometary stream modeling.
Simple analytic treatments produced some surprises as well.Nesvorný et al. (2011a) showed that implementing the nominal collisional lifetime of Grun et al. (1985) in their dynamical model resulted in the removal of the majority of radar-detectable meteoroids before they reach Earth; the authors estimated that collisional lifetimes 30 − 100× longer are necessary to reproduce the observed orbital distributions of meteors at Earth.Pokorný et al. (2014) implemented the Steel & Elford (1986) collisional model to the population of meteoroids originating from Halley-type comets, which is characterized by a broad range of orbital eccentricities.Pokorný et al. (2014) found a result like that of Nesvorný et al. (2011a): reproducing these long-period-comet meteoroids required collisional lifetimes 20 − 50× longer than those estimated in Steel & Elford (1986).Soja et al. (2019), in generating the latest version of the Interplanetary Meteoroid Environment Model, concluded that two different collisional lifetime categories provide the best fit to the Zodiacal Cloud measurements; for meteoroids with diameters D < 125 µm, the Grun et al. (1985) collisional lifetimes were found to be adequate, while for D > 125 µm, the authors required collisional lifetimes that are 50 times longer.Because smaller particles have shorter dynamical lifetimes (due to the greater importance of Poynting-Robertson drag at small sizes; Burns et al. 1979), they might be insensitive to collisions even for the nominal collisional lifetimes estimated from Grun et al. (1985).
The two most commonly used methods for estimating collisional lifetimes were developed in Grun et al. (1985) and Steel & Elford (1986).Both of these works assume a certain shape, density, and velocity distribution of particles in the Zodiacal Cloud based on observations available at the time, and calculate the collisional lifetime of a meteoroid on a given orbit.There are two major differences between these two approaches.First, Grun et al. (1985) characterized meteoroid strengths by a binding energy that the kinetic energy of a projectile must exceed in order to cause catastrophic destruction, whereas Steel & Elford (1986) only required a target-to-projectile radius ratio of ∼ 30 − 40.Second, the Steel & Elford (1986) collisional lifetimes are dependent on a meteoroid's orbital inclination, while in Grun et al. (1985) the collisional lifetimes are only a function of a meteoroid's heliocentric distance and mass.
Each approach has flaws.Disregarding the way in which impactor energy scales with heliocentric distance, as Steel & Elford (1986) do, can over/underestimate the rates of destructive collisions in the closest/farthest regions of the Zodiacal Cloud (Mann et al. 2004).Neglecting inclination in the manner of Grun et al. (1985) can result in underestimating the collisional lifetimes of meteoroids on highly-inclined orbits and overestimating the lifetimes of retrograde particles with orbits close to the ecliptic (Steel & Elford 1986).Note that Grun et al. (1985) and Steel & Elford (1986) provide similar estimates for collisional lifetimes for meteoroids with orbits close to the ecliptic at 1 au.Finally, any approach to collision modelling that is not fully three-dimensional underestimates the effects of mean motion resonances (Nesvold & Kuchner 2015b) and secular forcing (Nesvold & Kuchner 2015a) on collision rates.
To overcome the various drawbacks and assumptions of analytic methods to estimate the effects of mutual dust and meteoroid collisions, Stark & Kuchner (2009) devised a numerical algorithm for treating collisional evolution of any particle cloud.This new "collisional grooming" algorithm is a self-consistent method that uses a three-dimensional distribution of model particles and particle trajectories.The Stark & Kuchner (2009) algorithm does not make any assumptions about the shape of the dust cloud or the collisional lifetime but rather uses the modeled dust spatial and velocity distribution to estimate the rate of collisions at each point in the cloud.Using a collisionless dust cloud model as a starting point, the method iteratively applies a collision correction to the dust cloud density until a convergent steady-state solution is achieved: that is, the collision-less dust cloud groomed by the steady-state solution yields the same steady-state solution.The Stark & Kuchner (2009) method has been applied to model dust in the outer solar system (Kuchner & Stark 2010;Poppe 2016) where the quantity of model constraints is scarce.This method has not previously been applied to the inner, more easily observable, regions of the Zodiacal Cloud.
Motivated by the surprisingly large collision lifetimes required by previous models and by the abundance of different dust and meteoroid related constraints available for the inner solar system, we set ourselves the following goals: (1) create a publicly available version of the Stark & Kuchner (2009) code that can handle large dynamical dust models able to reproduce various inner solar system phenomena, (2) find the best fit to currently available constraints in the inner solar system using a new model of Jupiter-family comet (JFC) meteoroids, (3) compare the collisional lifetimes estimated from our new model to those used in previous studies, and (4) constrain the material strength of JFC meteoroids by using the collision rates of our best-fitting model to constrain their binding energy.

METEOROID MATERIAL CHARACTERISTICS AND COLLISION CONDITIONS
When a meteoroid (target) collides with a member of the particle cloud (projectile), this event has three possible outcomes: (1) the target suffers negligible damage; (2) the target incurs partial damage; and (3) the target is catastrophically destroyed leaving behind a cloud of daughter particles.In this article we only consider outcome (1) and (3), where we assume that the cloud of daughter particles is blown out of the solar system due to radiation pressure and thus can be neglected.We neglect the second scenario, the accumulation of damage on each particle, due to unknown effects of particle-particle collisions on the immediate structure of individual particles (as suggested in Grun et al. 1985).
In our analysis, we only check for catastrophic disruption as though the particle is a fresh particle, regardless of its dynamical age.The particle erosion effect would in general decrease the particles' collisional lifetimes, but due to the lack of any empirical/laboratory constraints we reserve that topic for future work.
Our catastrophic collision conditions follow the formalism from Krivov et al. (2005), in which a catastrophic collision occurs when the kinetic energy of the impactor E imp is greater than or equal to the target's binding energy E bind : where A s and B s are material constants, m tar is the mass of the target, R tar is the radius of the target, m imp is the mass of the impactor, v coll is the relative collision velocity, and Q * D is the critical specific energy.In this article, we use MKS units unless explicitly stated otherwise (e.g., for the semimajor axis we use astronomical units).Krivov et al. (2006) set A s so that Q * D = A s (R tar /1cm) Bs = 10 6 erg g −1 for a 1 meter target and set B s = −0.24.These choices correspond to A s = 3.02 × 10 6 erg g −1 in CGS units or A s = 100.0J kg −1 in MKS units.For comparison, Grun et al. (1985) uses the following expression for their binding energy: where S c is the unconfined compressive strength in kilobars (kbar), which is typically between 0.5-4 kbar or 50-400 MPa (Ostrowski & Bryson 2019).Comparing this expression to Eq. 2 assuming S 0.45

C
≈ 2, we see that Grun et al. (1985) used A s ≈ 300 J kg −1 , a value approximately three times larger than the Krivov et al. (2006) value (A s = 100 J kg −1 ).Rigley & Wyatt (2021) in their model for model for JFC meteoroid and dust cloud find values of A s = 660 and B s = −0.90.Their value of A s is comparable to those of Grun et al. (1985), however, their slope B s is much steeper as shown in comparison to other relevant works in their Figure 25.
After reorganizing a few terms in Eq. 2, we can obtain a ratio R max between the maximum destructable target radius R max,tar for a given impactor radius R imp and v coll : . (5) Figure 1 plots Eq. 5 for the range of impactor diameters modeled here and for collision velocities up to v coll = 50 km s −1 .

COLLISIONAL GROOMING CODE DESCRIPTION
Motivated by the findings of Stark & Kuchner (2009) and Kuchner & Stark (2010), we decided to implement a new version of the Stark & Kuchner (2009) collisional grooming algorithm that includes several optimizations and improvements.The entire code is written in C++ and is readily accessible at the project GitHub page1 together with a simple test case particle model described in this section.The code does not currently support multi-CPU/GPU processing.The basic collisional algorithm is described in Stark & Kuchner (2009), but we review it briefly here.

The Seed Model
As an input, the grooming algorithm requires records of the dynamical evolution of a particle cloud in Cartesian coordinates: for each particle, the code requires a list of the particle's position and velocity vectors as it evolves over its entire lifetime, from creation to destruction, assuming no collisions.Such a list is commonly obtained by simulating a population of particles in a numerical N -body integrator as we describe in Sections 4 and 5.1.
Each particle in an N -body simulation starts with some initial ⃗ r, ⃗ v and then dynamically evolves through a complex pathway until it meets its end in one of several possible particle sinks (e.g., gets too far from the central star, disintegrates when too close to the central star, hits a planet/moon, etc.).Our grooming code has no requirements for the end conditions; its only requirement is that the particle properties are recorded at uniform time intervals T record .Likewise, the code does not a priori expect that the particles' positions and velocities are recorded in either barycentric or asterocentric (star centered) coordinates, and can handle multi-star systems without any adjustments.This suite of particle records form the collision-less "seed" model.
Since the collisional grooming algorithm seeks a steady-state solution, the time domain and space domain merge.Each particle in the N -body simulation is replaced by a pathline representing many particles.Each of these pathlines is assigned a birth number density, n init , which decreases along its path due to collisions with other particles.
All the particle records are distributed into a three-dimensional grid of bins of equal dimensions, creating a three-dimensional histogram that represents the total particle density of the cloud.Each bin can be imagined as a box that contains recorded information of all particles that happen to be inside the box.Our code stores the particle's velocity vector ⃗ v, diameter D, and number density n in each bin, for all the particle pathlines that have records in that bin.

Grooming With Fewer Iterations
Once a collision-less seed model is loaded into the 3D histogram, we can investigate how collisions attenuate individual causal pathlines.When a particle with index i passes through a particular 3D histogram bin located at ⃗ r = (x, y, z) with the collisional depth τ coll (⃗ r), its pre-collision number density n pre,i is attenuated as where n post,i is the post-collision number density that particle with index i carries to the next record interval.We refer to the process of adjusting the number densities via Equation 6as "grooming".The value of τ coll effectively describes the number of collisions a particle experiences per sampling interval and depends on the model setup.
The collisional depth is calculated in each bin in our code as where n k is the particle density of the k th particle record in the bin, σ i,k = π(R i +R k ) 2 is the collisional cross section assuming that both particles have negligible escape velocities compared to their relative encounter velocities (no gravitational focusing), R i and R k are the radii of the investigated and k th particle respectively, v rel i,k is the relative velocity between the investigated and k th particle, where the sum is over all k records within an investigated bin.The collisional lifetime T coll -i.e., the collision e-folding lifetime -is defined as and describes the timescale on which the weight of the particle cloud is decreased to 1/e = 36.79% of the original weight; that is, 63.21% of particles are destroyed after one T coll .These approximations only work if the record sampling frequency is much greater than the average time between collisions, and if the bins are small enough to resolve any 3-D structures in the cloud, including those created by collisions.If the bins are too large, the code can over-estimate the rate of collisions and skew the final results.
All particle pathlines are processed independently: each one of them is individually groomed using the same initially collision-less particle cloud, which results in the first iteration of a collisionally processed particle cloud.However, since all particles are attenuated based on the collision-less model, this first iteration overestimates the collisional depths in each bin and arrives at a solution that underestimates the particle densities in each bin.Therefore, we perform further iterations of the collisional grooming process until the pre-and post-iteration state difference is smaller than a set precision.Stark & Kuchner (2009) showed that such an iterative process jumps between overestimated and underestimated states but can converge to a steady-state solution.
However, we find that this convergence can be slow or non-existent for dense particle clouds with high values of collisional depth τ coll .Therefore, we add an additional step to the process beyond the Stark & Kuchner (2009) algorithm.After each iteration, we calculate the difference measure D between the pre-and post-iteration number density for each particle in the particle cloud.Since we know the number density n of each particle before and after the collisional iteration step, we get where k iterates over all particle records in our 3D grid.To decrease the effect of over/underestimating the collisional depths of the cloud, we combine the pre-and post-iteration states of the particle clouds and assigning all particles in the post-iteration state the following number density: where Equation 11 expresses what weight the pre-and post-iteration state hold.We arrived at this prescription by testing various functional forms of (S) and various values of D. These different forms involved linear combinations of linear, quadratic, cubic, exponential, logarithmic, trigonometric, and hyperbolic functions.We find that the hyperbolic tangent used in Eq. 11 works best in various tests conducted in Sec. 4. Dust clouds with different particle distribution or densities could benefit from more customized re-weighting factors.If the final state of the collisionally groomed cloud is approximately known, then it can be prescribed in the first iteration to speed-up the computation significantly.
Compared to the original Stark & Kuchner (2009) setup, our approach converges 10-100 times quicker for thick clouds with frequent collisions (τ coll > 10 −3 ), when the difference for the cloud state before and after iteration is high, i.e. for D > 0. The algorithm also performs better than Stark & Kuchner (2009) even for particle clouds with fewer collisions; however, the improvement is less significant in these cases because they only tend to require only several (< 10) iterations.

SIMPLE TEST CASE
We tested our code on a small dataset, reproducing a test from Section 2.3.5 of Stark & Kuchner (2009).We simulated the evolution of a cloud of 10,000 particles orbiting a Sun-like star with no planets, where all particles are released in the form of a thin circular ring.For this particular seed model, the effects of collisions can be expressed analytically, as shown in Wyatt (1999).All particles in our model were D = 240 µm in diameter and had bulk densities of ρ = 2000 kg m −3 ; that is, the ratio of radiation pressure to gravity was β = 2.375 × 10 −3 .
We express the effects of radiation pressure and Poynting-Robertson drag on dust particles in terms of the parameter β, which is the ratio between the radiation pressure force F rad and the gravitational force of the central star F grav : where L ⋆ is the central star luminosity, L ⊙ = 3.828 × 10 26 W is the solar luminosity, M ⋆ is the mass of the central star, and M ⊙ = 1.98847 × 10 30 kg is the solar mass, Q pr is the radiation pressure coefficient that we assume to be unity (Q pr = 1).For a detailed description of the effects of stellar radiative forces on the dynamics of dust particles and meteoroids see Burns et al. (1979) or Murray & Dermott (1999).For the rest of this paper, we assume a density of ρ = 2,000 kg m −3 , the value used in Nesvorný et al. (2011a).Note that particle density influences both the particle-particle collision criteria as well as constrainable quantities such as the mass flux, particle cross-section, etc.
Using the same initial conditions as in Stark & Kuchner (2009), all particles in our test case model were released with an initial semimajor axis of a = 10 au, on circular orbits with eccentricity e = 0, and inclinations I = 10 • .Their initial longitudes of the ascending node Ω, arguments of pericenter ω, and mean anomaly M were picked randomly between 0 and 360 • .On their release, all particles are pushed to slightly more eccentric orbits due to radiation pressure.Then they spiral toward the central star via PR drag (Burns et al. 1979).We recorded the particle position and velocity vectors every 6956 years until the particles reached a heliocentric distance of 0.5 au.We used the SWIFT RMVS 3 numerical integrator (Levison & Duncan 2013) to create the seed model.For the collisional grooming, we selected a bin size of 0.05 au where bins extend up to 10 au from the central star.
Figure 2 shows the final distribution of normalized surface density Σ/Σ 10au as a function of heliocentric distance r hel produced by our new, improved collisional grooming algorithm.Our simple scenario closely resembles the one-dimensional analytic model in Eq. 3.36 of Wyatt (1999), which has a surface density of : where η 0 is the ratio between the PR decay time t PR and the collisional lifetime t coll at r 0 .For η 0 < 1, the collisional lifetime is longer than the PR drag timescale and the majority of particles spiral inward to the central star without being collisionally destroyed.On the other hand, for η 0 ≫ 1, collisions dominate the dynamical evolution of particles and only a small fraction of the original population survives long enough to reach small heliocentric distances.Figure 2 compares the analytic prediction from Eq. 13 (dashed black line) to the results produced by our collisional grooming code for the test case (gray solid line) for five different values of η 0 spanning 5 orders of magnitude.Figure 2 shows that our code correctly estimates the rates of collisions in this simple scenario and matches predicted profiles.
For η 0 = 1000.0,the particle cloud is very dense and the rate of collisions is approaching the resolution limit; you can see slight differences between the analytic model and the numerical model toward the center of the cloud.These errors can be alleviated by running the code on a finer grid with Figure 2. Normalized surface density Σ/Σ 10au of the steady-state solution as a function of heliocentric distance for our simple test scenario of dust migrating from heliocentric distance of 10 au via PR drag with no planets included (gray solid line).All collisions between particles are destructive in this case and the collision products are omitted.The analytic solution given by Eq. 13 is represented by a black dashed line.We ran our collisional grooming code for five different values of η 0 , which is linearly proportional to the collision-less cloud mass, surface area, or the total number of particles.The analytic solution agrees with the solutions of our code for all five values shown here.
more particles or with higher recording frequency (shorter T record ).Similar limitations were observed in Stark & Kuchner (2009).Such high values of η 0 can be found in massive debris disks such as the one orbiting β Pictoris, or for centimeter sized particles.However, in real-life scenarios the dust clouds are composed of particles with (A) a wide range of sizes and (B) different compositions, where both of these factors can change the value of η 0 by orders of magnitude even within the same dust population.

COLLISIONAL GROOMING OF JUPITER-FAMILY COMET METEOROIDS
The main goal of this article is to study how collisions affect the main component of the Zodiacal Cloud: the cloud of dust and meteoroids originating from Jupiter Family Comets (JFCs).Numerous independent models show that JFC meteoroids dominate the flux and number density of the Zodiacal Cloud over a wide range of particle sizes in the inner solar system (Koschny et al. 2019).Infra-red observations from IRAS and COBE were best reproduced when the major component of the Zodiacal Cloud beyond 1 au originates from JFCs (Nesvorný et al. 2010(Nesvorný et al. , 2011a)).Indirect observations of cosmic dust entering the Earth's exosphere point to a dominant JFC component as well (Carrillo-Sánchez et al. 2016, 2020).Multiple independent long-term observations of radar meteor orbit radars showed the dominance of a short-period cometary component in the sporadic meteoroid background observed at Earth (Campbell-Brown 2008;Galligan & Baggaley 2004;Janches et al. 2015).
In this section, we present a "seed" model of JFC meteoroids that is an extension of the Nesvorný et al. (2011a) model.We add collisions to this new dynamical model using our newly developed collisional grooming code.We compare this model to various independent observations to find a unique, best fit to the data.This process allows us to characterize the population of JFC meteoroids and the material parameters that govern their collisions, A s and B s (see Eq. 2).

Model of Jupiter-family comet meteoroids
To cover the large size range needed to reproduce our various constraints and the effects of collisions between dust particles and meteoroids originating in JFCs, we need to construct a new dynamical model.The models presented in Nesvorný et al. (2011a); Pokorný et al. (2018Pokorný et al. ( , 2019)); Soja et al. (2019) all suffer from several shortcomings: (1) they do not cover a wide range of particle diameters, where for collisional purposes we require particles as small as 1 micrometer and as large as 1 centimeter (10,000 micrometers); (2) they have disproportionate and sparse size bins; and (3) they select source particles to integrate that have been launched into bound orbits, which biases the model.
This pre-selection of particles on bound orbits is usually performed to obtain a better yield of numerical simulations with respect to the computational time.However, this improved yield comes at a cost.The blow-out heliocentric distance is r ⋆ = 2βa; thus, the particles removed from simulations by the preselection process tend to be smaller ones (Moorhead 2021).These small particles (β−meteoroids) can disrupt larger particles, so omitting them could severely distort the Size-Frequency Distribution (SFD).This practice might have caused earlier models to over-predict the number fluxes for micron-sized particles observed by various spacecraft (Malaspina et al. 2020).
To address shortcomings (1) and ( 2), we created a model that covers a large range of particle diameters D. We selected β values for our modeled population of JFC dust and meteoroids ranging from β = 0.837 (D = 0.6813 µm) to β = 5.7 × 10 −5 (D = 10, 000 µm) in 26 logarithmically spaced bins: To address shortcoming (3) of previous models (pre-selection into bound orbits), we launch particles from their parent bodies assuming that the parent bodies are uniformly spread in mean anomaly M , regardless of whether we expect the particles to remain bound.
Our model uses the same parent body distribution as Nesvorný et al. (2011a).Parent body orbits are assumed to be the same as those of objects originating in the Kuiper belt and transferred into the Jupiter-family comet orbits (Levison & Duncan 1997).The orbital elements are recorded/saved once they reach critical perihelion distance q ⋆ ; i.e., their perihelion distance q < q ⋆ .We use 10 different values of q ⋆ ∈ [0. 25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50] au.Cumulative distributions of semimajor axis a, inclination I, eccentricity e, and perihelion distance q of parent bodies in each q ⋆ bin are shown in Fig. 3.
Each meteoroid in our simulations is released from a randomly selected comet from the parent body population in the q ⋆ bin.The meteoroids are released with a, e, I corresponding to their parent bodies, while the remaining orbital elements Ω, ω, M are selected randomly between 0 and 2π.All meteoroids are assumed to be test particles and are subjected to the gravitational force of the sun and eight planets, the effects of radiation pressure, and PR drag.The integration time step is 1 day (86400 seconds), where close encounters with planets are handled separately by the SWIFT RMVS 3 integrator.The particles are numerically integrated until they are removed from the simulation once they 1) are close to the sun (r < 0.005) au; or 2) impact one of the eight planets; or 3) are too far from the sun (r > 10, 000) au.Meteoroids in our model can be released on initially hyperbolic orbits due to radiation pressure which corresponds to a more realistic meteoroid production scenario where the comets produce dust close to the Sun.We do not assume any dependence of particle production rate based on cometary orbit such as cometary activity models described in Jones (1995) or Crifo & Rodionov (1997) since Nesvorný et al. (2010) showed that cometary activity itself is insufficient to supply the current Zodiacal Cloud; cometary fragmentation, not solar-driven out-gassing, generates most of the dust (Rigley & Wyatt 2021).
In order to obtain the same number of particle records in each diameter range, we increased the number of simulated particles in our simulations until we recorded 10 GB for each particle size bin and 1 GB per perihelion distance bin.This amounts to 260 GB of data in binary format, which is approximately 8.2 × 10 9 records.Particles of each size and originating from different q ⋆ populations (Fig. 3) have vastly different dynamical timescales.This means that the initial number of particles in different size and perihelion bins can differ by orders of magnitude, since small particles disperse quickly from the system while the largest particles in our sample can evolve on timescales that are 5 orders of magnitude longer.Table 1 shows the total number of released particles per each size and perihelion distance cut-off bin.The amount of data provided by our JFC model is more than sufficient for the purpose of this work.We ran several tests and even 5% of the total amount of orbital records were sufficient to provide the same results as the full data set..6813 4,533,921 7,313,546 31,571,532 31,571,532 31,571,532 31,571,532 31,571,532 31,571,532 31,571,532 31,571,532 1.000 3,868,972 1,806,057 1,316,070 1,326,277 2,084,598 3,150,065 3,792,783 9,852,590 25,538,430 31,567,273 1.468 3,441,864

Observational Constraints and Model Parameters
Before we discuss the collisional grooming of our JFC dynamical model, let us introduce the observational constraints and the goodness-of-fit functions we will use to measure how well our model fits these constraints.The next five sections introduce the model constraints: 1. Heliocentric distance dust density profile (Sec.5.2.1) 2. Mass accretion rate at Earth (Sec.5.2.2) 3. Orbital distribution of meteors at Earth (Sec.5.2.3) 4. Total cross-section of the Zodiacal Cloud (Sec.5.2.4) 5. Size-frequency distribution at Earth (Sec.5.2.5)

Heliocentric distance density dust profile
The heliocentric distance (r hel ) profile observed by Helios space probes (Leinert et al. 1981) shows that the surface brightness of zodiacal light in the direction perpendicular to the eclipcic plane scales as S ∝ R −2.3±0.1 , which conversely means that the particle density scales with the heliocentric distance as n ∝ r −1.3±0.1 hel .These values hold for heliocentric distances from r hel = 0.3 au to r hel = 1.0 au and for z−distances close to the ecliptic: z < 0.05 au.These values were confirmed by the Parker Solar Probe (Howard et al. 2019).In our model, we calculate the dust cloud density n near the ecliptic by averaging the two bins adjacent to z = 0 for each heliocentric distance r hel in the model.Then we fit the single-power law to the density profile and obtain the model power law index M R and compare it to the value measured by Helios and PSP, O R = −1.3± 0.1.

Mass accretion rate at Earth
The mass of meteoroids accreted onto Earth has been estimated using various ground-based and space-borne experiments.Love & Brownlee (1993) estimated the accretion rate from impact rates on the Long Duration Exposure Facility (LDEF) satellite to Ṁ⊕ = 109, 500 ± 54, 800 kg day −1 .This value was later re-interpreted using the hydrocode iSALE in Cremonese et al. (2012) and the accretion rate was decreased to Ṁ⊕ = 20, 300 ± 2, 700 kg day −1 if the mass accretion is dominated by main belt asteroid meteoroids, and to Ṁ⊕ = 11, 500 ± 1, 400 kg day −1 if the mass accretion is dominated by cometary meteoroids.Moorhead et al. (2020) revisited the LDEF impact record and added the results from the Pegasus experiment that in 1965 measured the flux of dust particles on 194.5 m 2 of collecting area.Moorhead et al. (2020) were not able to reconcile both impact experiments, noting the significant uncertainties connected with both experiments and complexity of the correct interpretation of acquired data.
Other ways to constrain the meteoroid mass flux at Earth use metal abundances in Earth's mesosphere (Carrillo-Sánchez et al. 2016, 2020) or accumulation of ablated and non-ablated meteoric material in Antarctic ice-layers (Rojas et al. 2021).Like space-borne estimates, these methods are model-dependent.For instance, Carrillo-Sánchez et al. (2016) showed that using different sizefrequency distributions, reflecting observations from IRAS and Planck infra-red observations of the Zodiacal Cloud, provides ∼50% different estimates for the meteoroid mass flux at Earth.However, these approaches have the advantage that they can differentiate between different meteoroid populations impacting Earth.
We selected the mass flux of JFCs at Earth M M = 19, 600 ± 7, 500 kg day −1 from Carrillo-Sánchez et al. (2020) as our primary mass accretion constraint.This measurement provides the most recent estimate and agrees within uncertainty ranges with Rojas et al. (2021) who analyzed a large collection of melted and unmelted meteoroids from the Concordia site in Antarctica.Note that Carrillo-Sánchez et al. ( 2016) provides a different value for JFC mass flux at Earth, M Ṁ = 34, 600 ± 13, 400 kg day −1 , which we also consider in our analysis as an alternative value.

Orbital distribution of meteors at Earth
Two meteor orbit surveys have acquired more than 5,000,000 orbits: the Canadian Meteor Orbit Radar (CMOR; Campbell-Brown 2008) and the Southern Argentine Agile Meteor Radar (SAAMER; Janches et al. 2015).Additionally, several other facilities have acquired hundreds of thousands of meteor orbits, including the Advanced Meteor Orbit Radar(AMOR; Galligan & Baggaley 2004), and the Middle and Upper Atmosphere Radar (MU; Kero et al. 2012).Jones et al. (2001); Nesvorný et al. (2011a) showed that JFC meteoroids preferably feed the helion and anti-helion sources observed at Earth; i.e., the meteors appearing to come from sun-ward and anti-sun-ward directions.While other meteoroid populations are detected at Earth as well, they are prominently concentrated into the apex region (OCCs; Nesvorný et al. 2011b) and north/south toroidal region (HTCs; Pokorný et al. 2014); see also Jones (2004).The main belt asteroids release meteoroids that might appear in helion/antihelion sources.However, the impact directions of these meteoroids are scattered over a wide range of angles around the celestial sphere and difficult to detect by radars due to their relatively small Earth impact velocities (Wiegert et al. 2009).The main belt asteroids dominate the number flux of sporadic centimeter+ sized meteoroids (D > 10, 000 µm) at Earth (Shober et al. 2021).However, that size range is currently beyond the scope of our model; we can assume that the meteors observed by meteor radars in helion/anti-helion sources are dominantly sourced from JFCs.
One of the significant differences between CMOR and SAAMER is in their transmitted power and the beam pattern, which translates into a difference in minimum detectable meteoroid mass.SAAMER, with its 60 kW transmitter and more focused beam, can see approximately an order of magnitude smaller particles in diameter than CMOR's 6/12 kW antennas with a wide beam pattern (Weryk & Brown 2012).Smaller particles suffer fewer collisions during their dynamical lifetimes than their larger counterparts (see e.g., Pokorný et al. 2018, for different meteoroid sizes at Mercury).This effectively means that while SAAMER provides valuable data about meteoroid environment, it cannot constrain the effects of collisional grooming as efficiently as CMOR.Nesvorný et al. (2011a) showed this contrast between AMOR and CMOR.AMOR with its narrow beam had a similar sensitivity to SAAMER and recorded meteors that showed almost no signs of collisional erosion.AMOR meteors were recorded with low eccentricity orbits resulting from efficient PR drag.On the other hand, CMOR meteors were mostly recorded having highly eccentric orbits.This shows that larger meteoroids detected by CMOR are collisionally removed before their orbits are circularized via PR drag and only the dynamically young population with high-e orbits is detected at Earth.
For these reasons we use Campbell- Brown (2008) helion and anti-helion raw orbital element distributions as our primary meteor orbit distribution constrains.These distributions were collected from approximately 1,200,000 meteors observed at these two sources.The CMOR helion source is defined in the heliocentric ecliptic longitude λ−λ ⊙ and latitude β as: The sensitivity of meteor radars is complex in nature (e.g.Janches et al. 2015) and involves multifacility observations, extensive analysis, and numerous assumptions (Weryk & Brown 2013).For the purpose of this paper, we adopt the ionization limit used in Campbell- Brown (2008) as where the mass of the meteoroid m met , the relative hyperbolic excess impact velocity v ∞ , the escape velocity v esc , and the impact (gravitationally focused) velocity v GF are in SI (MKS) units, i.e, a meteoroid with m met = 10 −7 kg and v GF = 30, 000 m s −1 has I ⋆ = 1.Campbell- Brown (2008) estimated that I ⋆ = 1.0 is the detection limit for CMOR, which is the value we use in this paper.Note that Campbell- Brown (2008) derived their expression using an extensive review done in Bronshten (1983) and our Eq.15 provides a simplified description of meteor ionization.The process depends on meteor velocity, composition, or atmospheric conditions (Jones 1997) and is still an open topic of research (e.g., Brown & Weryk 2020).
In our model, we first calculate the orbital elements of each meteoroid a, e, I using the heliocentric position and velocity vectors from the numerical simulation rather than the binned data stored in the 3D histogram.Then we calculate the collisional probability C p of each meteoroid record with an idealized Earth (a = 1.0 au, e = 0.0, I = 0.0 • ) using the Kessler (1981) method.This method uses the orbital elements of the projectile and the target to analytically determine their collisional probability per unit time.For our idealized Earth this method is identical to that of (Opik 1951).Multiplying C p by the size dependent model normalization parameter (Sec.5.2.5 yields the number of impactors at Earth.Using the particle size, SFD index, v GF and Eq. 15, we determine what fraction of meteoroids in the size bin impacting Earth is detectable by the radar; i.e., I ⋆ ≥ 1.After investigating all meteoroid records in the numerical simulation, this procedure results in four histograms of a, e, I, and v ∞ for meteoroids that would be detectable by CMOR, which we collectively summarize as our model measurement M O .The observed distribution of orbital elements from CMOR radar is then O O .5.2.4.The total cross-section of the Zodiacal Cloud Gaidos (1999) derived an effective emitting area for the Zodiacal Cloud of Σ ZC = 1.12 × 10 17 m 2 by assuming that the dust particles emit blackbody radiation at 260 K (Reach et al. 1996) and have the bolometric luminosity is 8 × 10 −8 L ⊙ = 3.12 × 10 19 W (Good et al. 1986).Nesvorný et al. (2010) used a dynamical model constrained by Infrared Astronomical Satellite (IRAS) observations of thermal emission to estimate Σ ZC = 2.0 ± 0.5 × 10 17 m 2 .Nesvorný et al. (2011a) estimated the total crosssection of their modeled Zodiacal Cloud to be 1.7 × 10 17 < Σ ZC < 3.4 × 10 17 m 2 , in agreement with their previous estimate, but approximately twice the Gaidos estimate.Unable to decide between these two estimates, we decided that for the purpose of this paper we will conservatively assume two scenarios: (A) Σ ZC = 2.0 × 10 17 m 2 , based on dynamical modeling works, and (B) Σ ZC = 1.12 × 10 17 m 2 using Gaidos (1999).In our model, we calculate the total cross-section of the Zodiacal Cloud by summing the cross-sections of all meteoroids in the model assuming they are perfect spheres; i.e, M Σ = i 0.25πD 2 i .

Size-frequency distribution at Earth
The interplanetary dust particles and sporadic meteoroids considered in our model range in size from D = 1 µm up to D = 10, 000 µm.We exclude particles smaller than 1 µm as a significant fraction of sub-micron particles that intercept the Earth are interstellar (Strub et al. 2019).We place our upper limit at 10, 000 µm, or 1 cm, well below the lower size limit for asteroids (∼1 m) and approximately the size above which meteoroids are too scarce to be detected in large numbers by ground-based networks.Even within these limits, no single instrument can characterize the entire SFD of interplanetary dust particles and meteoroids.
Fortunately, at least for some range of meteoroid diameters the data is abundant enough to constrain the SFD at Earth.Pokorný & Brown (2016) determined the SFD of radar and optical meteors using more than 5,000,000 meteor echoes observed by CMOR and 3,106 optical meteors observed by CAMO (the Canadian Automated Meteor Observatory).Both samples were scrubbed of meteor showers to avoid contamination from meteoroid populations originating in singular sources (for daily fluctuations of the SFD see Fig. 7 in Pokorný & Brown 2016).From these two independent observations, the authors determined that the SFD for meteors with masses 10 −8 kg < m met < 10 −4 kg has a differential mass index of s = −2.1 ± 0.08.These values can be converted to a differential size-frequency index S = −4.3± 0.24 for the range of diameters 212.2 µm < D < 4571 µm assuming spherical particles and a bulk density of ρ = 2,000 kg m −3 .Galligan & Baggaley (2004), using the AMOR radar, showed that S is shallower for smaller meteors and reported S = −4.09± 0.03 for meteors as small as m met = 10 −10 kg (or D = 45.7 µm).
In our model, we fit the size-frequency distribution at the Earth with a broken power law defined by two differential size indexes a 1 and a 2 , the break point µ, and the offset b as: where 10 b is a normalization constant.The index a 1 defines the slope of the broken power law for particles with D larger than the break point µ, whereas a 2 represents the slope for the smaller particles in the SFD.As break points µ > 212.2 µm have not been reported, we will constrain the value of µ to be less than 212.2 µm.The value of a 1 is the observed differential SFD index for large particles -M α = a 1 -which we will compare to the observed value O α = −4.3± −0.24 from Pokorný & Brown (2016).
We are not aware of any recent observations that constrain the value of a 2 .Fixsen & Dwek (2002) derived a value of a 2 = −2.3± 0.1 and a SFD break point µ = 30 µm from DIRBE and FIRAS instruments on the COBE satellite that observed the spectrum of the Zodiacal Cloud.This value is not directly comparable with the SFD of meteoroids impacting Earth due to the different effect of gravitational focusing on smaller and larger particles.Therefore, we will only constrain the value of a 1 and discuss the SFD of particles with D < µ in our discussion section.

Goodness-of-fit function for model constraints
To quantify which of the model's free parameters provide the best representation of observed reality, we construct a goodness-of-fit function χ 2 .The combination of free parameters that minimizes χ 2 will be considered our best fit.We define the total model χ 2 tot as follows where χ 2 R , χ 2 Ṁ , χ 2 O , and χ 2 α are the goodness-of-fit of the heliocentric distance profile, mass accretion rate at Earth, orbital distributions of meteors at Earth, and size-frequency distribution index.These goodness-of-fit functions are defined as follows: where M i are the model values, O i are the observed values, and σ i are the uncertainties associated with the observed values for i ∈ {R, Ṁ , Σ, α}.Unfortunately, the orbital element distributions from Campbell- Brown (2008) do not come with uncertainty estimates.To calculate χ 2 O , we adopt from Campbell- Brown (2008) that from the total number of observed meteors N CMOR = 2.35 × 10 6 approximately P HE/AH = 51% are in the helion/anti-helion source.Campbell- Brown (2008) published distributions of four meteor orbital parameters: semimajor axis a, eccentricity e, inclination I, and the hyperbolic excess impact velocity v ∞ , each distributed into 100 uniformly spaced bins.All distributions in Campbell- Brown (2008) are normalized, so we estimate the average variance in each bin as Then we can estimate χ 2 O using following expression where the factor 400 comes from the summation over 400 histogram bins.Note that we average over the 400 orbital element bins to assess the average chi-square for the distribution.Otherwise, χ 2 O would overwhelmingly dominate χ 2 tot and the constraints offered by other observations would essentially be ignored.There are other sources of uncertainties connected with meteor radar measurements such as the ionization efficiency variability, uncertainties in meteor velocity, and direction determination (e.g., Mazur et al. 2020) that are currently beyond the scope of the paper.We ran simple tests by increasing and decreasing our σ 2 O by a factor of 2 and found no significant changes in the results.

Model parameters and minimization process
Loading the JFC meteoroid model into our collisional grooming code distributes 371,247,921 unique records into 26 logarithmically scaled size bins.We perform a search for the minimum value of χ 2 tot (defined in Eq. 17) over a four-parameter phase space (P α , P γ , P As , P Bs ) composed of: (1) the differential power law index for the SFD of the source population of JFCs meteoroids, P α ; (2) the JFC initial perihelion weighting parameter P γ , where each meteoroid stream carries weight W = q ⋆Pγ ; (3) a parameter that determines the material property constant A s , P As = 100A s J kg −1 (where A s is the parameter in Krivov et al. 2005, after converting cgs to MKS; see Section 2); and (4) the material property constant B s = P Bs (see Eq. 2).Parameters P α and P γ control the initial size-frequency and heliocentric distance distribution in the particle cloud, whereas P As and P Bs control the collisional grooming process.We use a factor of 100 in P As = 100A s to directly see the multiplication factor against values of A s used in Krivov et al. (2005) and Kuchner & Stark (2010).For additional discussion of the influence of P γ see Section 2.1 in Nesvorný et al. (2011a).
The absolute scaling of the number of dust particles in the Zodiacal Cloud is determined by the initial value of the total Zodiacal Cloud cross-section Σ ZC, init for the collision-less cloud.In Scenario A, we select Σ ZC, init = 2.0 × 10 17 m 2 , while for Scenario B we set Σ ZC, init = 1.12 × 10 17 m 2 .As soon as the collisions are iteratively taken into account Σ ZC decreases.Once a steady state is reached, the model has a total Zodiacal Cloud cross-section of Σ ZC, steady .The entire model is then re-scaled by a factor of Σ ZC, init /Σ ZC, steady > 1.0, and the model is collisionally groomed again until a second, augmented steady state is achieved.If the second augmented steady state is within 10% of Σ ZC, initi.e., | ln(Σ ZC, init /Σ ZC, steady )| < 0.0953 -then the collisional grooming is considered complete.If this condition is not satisfied, further rescalings are performed.A single rescaling is usually sufficient due to the tenuous nature of the current Zodiacal Cloud.

Additional constraints and observations
To keep this work concise, we include observational constrains that we are the most familiar with.Additional constraints for the shape of the Zodiacal Cloud and its Jupiter-family comet component are for example: the shape and brightness of the inner Zodiacal Cloud (Hahn et al. 2002); the Zodiacal Cloud emission spectrum and the size-frequency distribution break point at 14 and 32 µm (Fixsen & Dwek 2002); micro-cratering records on the Moon (Hoerz et al. 1975); the Parker Solar Probe dust count record from 0.05 au to 1 au (Malaspina et al. 2020); or the spacecraft impact crater records from Pegasus and LDEF (Moorhead et al. 2020).We compare our model results with the Parker Solar Probe in Sec.7.1 and leave the inclusion of more observational constrains for future work, where we will include additional major dust and meteoroid populations.

RESULTS
In this section, we insert the entire JFC dynamical model into the collisional grooming code and look for the model realization that best fits reality as determined by the constraints described in Sec.5.2.

Constraining the fitting parameter space -Scenario A
Due to the computational demand of each model realization (6-8 hours of real time) and the broad parameter space, we first focused on determining the influence of each fitting parameter separately.We fix three of the four different fitting parameters P in our initial reconnaissance and let the fourth vary across a wide range of values.After experimenting with 256 different permutations of four free parameters we found a solution that provided a satisfactory total goodness-of-fit χ 2 tot as well as partial goodness-of-fit values (χ 2 R , χ 2 Ṁ , χ 2 O , χ 2 α ) for P α = 4.47, P γ = 0.15, P As = 5500, and P Bs = −0.24.For the final round of our reconnaissance, we explore the effect of each individual parameter P while keeping the other parameters fixed at the previously mentioned values.Figure 4 shows variations in our four goodness-of-fit function values as a function of the four parameter values.Variations of χ 2 values of individual parameters allow us to define a focused fitting parameter space (gray areas in Figure 4).This parameter space constrains the ranges of all model parameters P and is used to fit all four model parameters at the same time.

Influence of P α
First, we examine the influence of P α , which represents the size-frequency distribution of particles at the source.Figure 4A shows the most significant drop in χ 2 tot for the differential size index at Earth (χ 2 α ) and the JFC mass accreted at Earth (χ 2 Ṁ ).These two functions go from indicating very poor fits (χ 2 > 10) to indicating very good fits (χ 2 < 0.1) inside the range we investigated (P α ∈ [3.6, 5.3]).The minimum values for these two goodness-of-fit functions are located between P α = 4.20 (for χ 2 α ) and P α = 4.50 (for χ 2 Ṁ ).We also see a significant drop in χ 2 R but only in the log-scale since the value is very close to 0. The value itself does not go above χ 2 R > 0.2 because we kept the value P γ = 0.15 which provides a good fit to the heliocentric distance distribution of JFC meteoroids in our model (Fig. 4C).The meteor orbit distribution is quite insensitive to P α as the radar detects meteors within a certain size range and thus, if the SFD is anchored at that size, changing the slope of the size-frequency distribution has a minimal effect on the distribution of detected meteors.Based on our initial analysis we select P α ∈ [4.2, 4.6] for our further analysis (gray region in Fig. 4A).

Influence of P As
Parameter P As strongly influences the collisional lifetime of each particle by making particles more or less susceptible to collisions.Since one of our model constraints is the total cross-section of the currently observed Zodiacal Cloud (Section 5.2.4), we can expect that P As mainly influences the total mass of particles impacting Earth, their orbital element distribution, and their size-frequency distribution.We do not expect this parameter to alter the heliocentric particle density profile significantly, because this profile is dominated by smaller particles that are evolving quicker than their larger counterparts.Figure 4B shows our initial prediction, where we see that for P As < 10 the goodness-of-fit functions are well-above χ 2 > 1 but improve significantly as P As increases.Once P As reaches a value of 1000, the goodness-of-fit values start to plateau, showing that collisions in the particle cloud are less and less important above this threshold.Once reaching P As > 10, 000, the collisions in the cloud rarely serve to destroy particles, and our constraints cannot distinguish between different models.All models are equivalent to a collision-less model above this threshold.This means that for P As > 10, 000, particles in our model are not sensitive to collisions.Based on our initial analysis we select P As ∈ [1000, 8000] for our further analysis (gray region in 4B).

Influence of P γ
As suggested by its definition, the parameter P γ most significantly influences the heliocentric distance profile of JFCs in the inner Zodiacal Cloud (Fig. 4C).Furthermore, this parameter changes the weighting of low-eccentricity and high-eccentricity orbits of meteoroids launched from JFCs because particles emitted near perihelion are ejected into more eccentric orbits.As a result, the eccentricity distribution of those meteoroids that impact Earth will vary with P γ .Despite this, Figure 4C shows that P γ has very little influence on χ 2 O and χ 2 Ṁ .This is due to the fact that χ 2 O and χ 2 M are quantities that are more sensitive to circularized meteoroids at 1 au and are accentuated by gravitational focusing at Earth.However, P γ does, as expected, strongly influence the value of χ 2 R , which has a minimum around P γ = 0.15.P γ also influences the SFD of meteoroids impacting Earth and thus the value of χ 2 α .While the minimum of χ 2 α lies around P γ < −0.50 for the default combination of parameters (P α = 4.47, P As = 5500, P Bs = −0.24),we previously showed that χ 2 α also strongly depends on P α .We observe that increasing P γ always steepens the SFD of the population of Earth's impactors.Since P γ is the only parameter that strongly changes χ 2 R we selected the interval P γ ∈ [0.125, 0.175] for our further investigation.

Influence of P Bs
Parameter P Bs controls the scaling of the meteoroid binding energy E bind with its radius (Eq.2).E bind is linearly proportional to A s and the meteoroid's mass, whereas B s determines the relationship between the binding energy and the radius of the target particle (E bind ∝ R Bs tar ).For B s < 0 and meteoroid radii < 1 meter (i.e., all meteoroids considered here), meteoroids are more difficult to break when compared to the scenario in which B s = 0, as their specific (per-unit-mass) binding energy increases.The fact that meteoroids are more resilient to fragmentation with decreasing size is something we would expect based on the latest in-situ observations of dust building blocks from ROSETTA (Güttler et al. 2019).A suite of instruments on-board the ROSETTA spacecraft observed dust particles and meteoroids of various sizes and morphologies and concluded that agglomerate dust particles are built-up from smaller sub-structures, with the smallest building block sizes commonly around 0.5 − 2.0 µm.Apart from completely rock-solid dust grains, the aggregate dust particles tend to be more loosely packed together with increasing size and the smaller number of sub-structure connections decreases overall meteoroid strength.
Based on our initial 256 parameter permutation analysis (Sec.6.1) and the targeted follow-up analysis, we find that P Bs has little influence on the goodness-of-fit functions below P Bs < −0.10 (see Fig. 4D).For P Bs > 0.0, the size-frequency distribution started to show waviness, an oscillation of the SFD index between D = 10 µm and D = 1000 µm.This results in a SFD that is impossible to fit with a broken power law and invalid values of χ 2 α for P Bs > 0.2.From this analysis, we conclude that B s does not have a strong effect on our goodness-of-fit as long as P Bs < −0.1.Therefore we decided to use the original Krivov et al. (2006) value, B s = −0.24,for the rest of the models in this article.

Finding the best fit
We ran 216 different model realizations within our fitting parameters constraints; the results are summarized in Figure 5A.In order to display four-dimensional information in an image, we represent the parameter P γ using three different symbols.Results of our focused χ 2 tot minimization in Scenario A follow the results of our one-dimensional analysis described in Section 6.1.Parameter P As , which scales the particle binding energy, produces small variations in χ 2 tot for P As > 2000.The initial JFC differential size index P α minimizes χ 2 tot around 4.35 ≤ P α ≤ 4.4 (denoted by the black rectangle).The minimum value of χ 2 tot = 1.13 was recorded for (P α ,P As ,P Bs , P γ ) = (4.40,8000, −0.24, 0.125), where all values in the black rectangle have χ 2 tot < 1.20 and thus within 6% of the minimum value.Figure 5B provides a closer look at the best fitting parameter combinations in our analysis.We constrained the color-scale to values with log 10 χ 2 tot < 0.1; i.e., χ 2 tot < 1.26.We see that the bestfitting models have P As > 6000 and P α = 4.4, where P γ does not have a significant influence on the results.
Our analysis shows that χ 2 tot > 1 for all parameters in our multi-dimensional parameter search, while in the one-dimensional analysis we were able to find parameter combinations for all four goodnessof-fit functions with values χ 2 < 0.1, except for the radar orbit distribution χ 2 O that we show is only sensitive to P As .The reason for the low-quality fit is that the individual goodness-of-fit functions for mass accreted at Earth χ 2 Ṁ and the SFD slope χ 2 α do not show their minimum values for the same combinations of parameters (Fig. 5C-D).This suggests a potential conflict between our constraints that our JFC model might not be able to reproduce for any combination of model parameters.For P α = 4.2, where we obtain the best fit for χ 2 α , the mass accreted at Earth is too high with Ṁ ≈ 34, 000 kg day −1 , whereas for P α = 4.45, where χ 2 Ṁ is the lowest, the differential index α is too high with a 1 ≈ −4.56.Now, we take a closer look at the ability of the best fit model to reproduce our constraints.Figure 6A-D shows that our model fits the radar meteor orbit distribution from CMOR (orange solid lines) well with the exception of a lack of meteors faster than v ∞ = 40 km s −1 .The remaining orbital elements observed in helion and anti-helion sources are very close to model predictions (black solid lines).The observation-model residuals (blue solid line) show expected noise due to the radar orbit uncertainty that washes out the fine structures seen in our model.For example the radar data are not able to capture the effect of inner mean-motion resonances with Jupiter while the model shows a dip at 2.5 au corresponding to the 3:1 mean-motion resonance.The lack of faster meteors can be attributed to long-period comet populations such as the Halley-type comets not included in our model (Pokorný et al. 2014) that are a part of the helion/anti-helion complex observed by CMOR.
The model cumulative size-frequency distribution shown in Fig. 6E (black solid line) spans over four orders of magnitude in particle diameter (or 12 orders of magnitude in particle mass).We fit the four-parameter broken power law (Eq.16) to this distribution (orange dashed line in Fig. 6E) where the fitting parameter a 1 represents the cumulative size index for the larger portion of particles.The observed differential size-frequency index O α = −4.3± 0.24 is shallower than our best model predicts M α = a 1 = −4.505,which is expected from our χ 2 analysis shown in Fig. 5.To match the observed values we would need the size-frequency distribution at source to be around a 1 = −4.2(i.e.P α = 4.2).For comparison we also show the lunar (blue dashed line) and interplanetary (pink dashed line) cumulative particle flux as estimated in Grun et al. (1985), where we multiplied the Grun et al. (1985) estimates by a factor of 3 to take into account the effect of gravitational focusing that JFC meteoroids experience at Earth (Pokorný et al. 2018).The gravitationally focused Grun et al. (1985) flux agrees well with our model between 80 < D < 10, 000 µm, where for smaller diameters our model predicts 1-2 orders of magnitude more particles than the Grun et al. (1985) model.This is due to our assumption of a single power law at the source, while there could be a break in the power law.We discuss the potential shortcomings of our single power law distribution in the Discussion section.
The variation in particle density with heliocentric distance close to the ecliptic (black solid line) shows a very good match to the observed power law scaling (orange dashed line; Fig. 6F).We fit the power law to our model values between 0.3 and 1.0 au.The model shows a decrease of the power law slope below r hel = 0.1 au due to the more frequent collisions in the innermost parts of the solar system.The rapid decrease in the particle density beyond 5.0 au is due to the initial distribution of JFCs.The particle density beyond 5.0 au is not representative of the overall dust distribution due to contributions of outer solar dust populations originating from long-period comets and Edgeworth-Kuiper belt objects (Poppe 2016).Note that the particle density in the innermost parts of the solar system (r hel < 0.1 au) might not be correctly represented due to currently unknown effects on dust that are being reported by the Parker Solar Probe (Stenborg et al. 2021) and that are not included in the model.
In Scenario A, we use the normalization to the Nesvorný et al. (2010) and Nesvorný et al. (2011a) model values, Σ ZC = 2.0 × 10 17 m 2 for the total cross-section of JFCs.Our results show that our model can either fit the observed mass accretion at Earth or the observed size-frequency distribution of meteoroids at Earth, but it fails to fit both constrains simultaneously.The radar meteor orbit distribution and the particle density variations with heliocentric distance are reproduced well for a variety of model parameters due to their insensitivity to parameter P α which determines the sizefrequency distribution at the source.From our analysis, we see that in order to fit all constrains well we would need to decrease the mass accreted at Earth by a factor of 2. In Section 5.2.4 we consider two values for Σ ZC : Scenario A: Σ ZC = 2.0 ± 0.5 × 10 17 m 2 ; Scenario B: Σ ZC = 1.12 × 10 17 m 2 , where their ratio is 1.79; a value very close to a factor of 2. We could also challenge the Carrillo-Sánchez et al. ( 2020) JFC mass flux and use the value estimated in Carrillo-Sánchez et al. (2016) who estimated O Ṁ = 34600 ± 13800 kg day −1 .Such a modified value allows for a very good fit at α = −4.2 and thus reconciling the inability to fit both the size-frequency distribution and the mass accretion rate at Earth at the same time.Using O Ṁ = 34600 ± 13800 kg day −1 we obtain the best fit for: P α = 4.20, P γ = 0.175, P As = 8000, P Bs = −0.24,with the goodness-of-fit value χ 2 tot = 0.298, indicating a very good fit to all constraints.Motivated by results from Scenario A, we now change the model total Zodiacal Cloud cross-section estimate to that of Gaidos (1999), i.e.Σ ZC = 1.12 × 10 17 m 2 .This will effectively decrease the amount of particles in our model to 56% of our first model in Scenario A and therefore decrease the collision rates between meteoroids by the same amount.This change should result in a decreased mass flux at Earth for the same model parameters and longer collisional lifetimes of all particles.

Scenario B: Constraining the fitting parameter space
We perform the same initial parameter space analysis as in Scenario A by fixing three of our fitting parameters and exploring the remaining fitting parameter for a wide range of values (Fig. 7).In order to reflect lessons learned in Scenario A, the fixed parameter values are the following: P α = 4.22, P γ = 0.15, P As = 5500, P Bs = −0.24.A significant difference compared to Scenario A is a much better agreement of minimum values of χ 2 Ṁ and χ 2 α for similar values of parameter P α in Fig. 7A.Both of these goodness-of-fit functions show a global minimum around P α = 4.2, showing a good sign that the reduced amount of particles in the cloud can provide a much better fit to observed constraints.As expected, the good fits are found for smaller values of parameter P As which does not show significant improvement beyond P As > 3000 (Fig. 7B).Interestingly, χ 2 R has a minimum for slightly larger value of P γ than in Scenario A, which is correlated with the shift of the fixed parameter to P α = 4.22 from P α = 4.47 in Scenario A (Fig. 7C).Parameter P Bs shows a worsening for all χ 2 functions for P Bs > −0.10, thus we keep it at P Bs = −0.24using the same value as in Scenario A.

Scenario B: The best model fit
For Scenario B, we ran 360 different model realizations using the parameter space defined in Sec.6.2. Figure 8 shows the results of our search.Compared to Scenario A, the total goodness-of-fit is much better with the minimum χ 2 tot = 0.292 for (P α ,P As ,P Bs , P γ ) = (4.20,5000, −0.24, 0.20) than that in our Scenario A analysis (χ 2 tot = 1.13).This is due to the fact that both the size-frequency distribution fit χ 2 α and the mass accreted at Earth fit χ 2 Ṁ have minima for the same location in our four-dimensional parameter space (Fig 8C -D).For P α = 4.20 and P As > 3000 both goodness-of-fit functions (χ 2 α ,χ 2 Ṁ ) are < 0.015, indicating a very good fit to the observed constraints.Since both χ 2 α and χ 2 Ṁ are very close to zero, the value of χ 2 tot is dominated by the goodness-of-fit to radar meteor orbital distributions with χ 2 O > 0.27 for all parameter permutations in our model.Note, that the quality of our best fit is similar to that found in Scenario A using the Carrillo-Sánchez et al. ( 2016) value for the JFC mass flux at Earth.Both of these best fits are achieved by using almost identical model parameters, where P α = 4.20, P As > 3000, 0.175 < P γ < 0.20.
Figure 9 shows the ability of our best fitting model to reproduce our observational constraints.From individual χ 2 values, we see that the Scenario B model provides a better fit in all analyzed aspects compared to the best fitting model in Scenario A. While the fit to meteor radar data and the particle density variations with the heliocentric distance is comparable to values recorded in Scenario A, the fits to the mass accreted at Earth and the size-frequency distribution index are reproduced very well with χ 2 ∼ 0.01.The shallower size index also provides a better match to Grun et al. (1985) for micron-sized meteoroids (in Scenario A the number flux at Earth was 2-3 orders of magnitude higher in our model compared to that in Grun et al. (1985) for D < 5 µm).

Summary of both fitting Scenarios
To summarize our model minimization analysis, we find that Scenario A provides a low quality fit χ 2 tot = 1.13 assuming the Carrillo-Sánchez et al. (2020) JFC mass flux at Earth.The model is able to either reproduce the size-frequency distribution at Earth for P α = 4.2 but provides ∼2× larger mass flux of JFC meteoroids at Earth, or reproduces well the mass flux for P α = 4.45 but then provides a too steep meteoroid SFD at Earth.This discrepancy can be reconciled by using the Carrillo-Sánchez et al. ( 2016) JFC mass flux estimate that provides a much better fit to all constraints and results in the total goodness-of-fit value χ 2 tot = 0.299 for the parameter combination (P α ,P As ,P Bs , P γ ) = (4.20,8000, −0.24, 0.175).
In Scenario B, we assume the total cross-section of meteoroids is 44% smaller than in Scenario A following the Gaidos (1999) estimate.In this case, we find that our model can fit all constraints presented here with χ 2 tot = 0.292 for the parameter combination (P α ,P As ,P Bs , P γ ) = (4.20,5000, −0.24, 0.20).Should we assume that the Carrillo-Sánchez et al. ( 2016) JFC mass flux is the correct value, Scenario B suffers from similar problems we encountered in Scenario A.
From our parameter search, we conclude that our best model fit is for the following parameters: P α = 4.20 ± 0.10, P As = 5000 ± 4000, P Bs = −0.24,and P γ = 0.20 ± 0.025.We estimate the 1 − σ interval for each parameter as the range where all individual χ 2 values are below unity.This means that Jupiter-family comet meteoroids are released from their source population with an initial differential size-frequency index α = −4.20 ± 0.10, their binding energy is E bind = 5 × 10 5 ± 4 × 10 5 R −0.24 met m met , and their weighting follows the trend q ⋆0.20 .

COLLISIONAL LIFETIMES
Having found the best fit to our constraints in Scenarios A and B, we can now calculate the collisional lifetimes for particles of various sizes and orbits in our model and compare them to other estimates used in recent Zodiacal Cloud models.In this Section, we compare the collisional lifetime T coll (Eq.8) obtained from our Zodiacal Cloud model to estimates from Steel & Elford (1986) and Grun et al. (1985): T coll (SE86) and T coll (G85).Each of these estimates has been used in numerous dynamical modeling papers.T coll (SE86) was used for example in Wiegert et al. (2009), Pokorný et al. (2014), Pokorný et al. (2018), Pokorný et al. (2019), and Yang & Ishiguro (2018) while T coll (G85) was implemented in Nesvorný et al. (2010), Nesvorný et al. (2011a), Soja et al. (2019), and Yang & Ishiguro (2018).
Figure 10 shows collisional lifetimes calculated in our best fit model in Scenario B using Eq. 8 for a meteoroid with diameter D = 50 µm assuming three different semimajor axes: a = 0.3 au (Fig. 10A), a = 1.0 au (Fig. 10B), and a = 3.0 au (Fig. 10C).The collisional lifetime was calculated by averaging the collisional lifetimes over the full range of the meteoroid's argument of pericenter ω, the longitude of the ascending node Ω, and the mean anomaly M , where for each of the three variables we took 100 values in the range (0, 2π), i.e. using an average of 10 6 different points in the phase-space (ω, Ω, M ).We calculated the collisional lifetime in our model for the full range of possible bound orbits, i.e. eccentricity 0 < e < 1, and inclination 0 < I < 180 • .
The collisional lifetimes in our model range over 5 orders of magnitude between the longest-lived particles and the shortest-lived ones, which are at retrograde and nearly parabolic orbits.Meteoroids on low-e and low-I orbits have the longest T coll and the collisional lifetime decreases with increasing e and I.This behavior appears at all semimajor axes.
Perhaps these variations in collisional lifetime are to be expected.The circular coplanar meteoroids orbit together with the bulk of the Zodiacal Cloud and therefore experience the smallest rates of destructive collisions due to very small relative impact velocities.On the other hand, retrograde meteoroids orbit against the flow and have much higher relative impact velocities.With increasing eccentricity, meteoroids penetrate into denser regions of the Zodiacal Cloud which naturally decreases their collisional lifetimes via increased rates of collisions.Moreover, their relative impact velocities are additionally increased due to higher orbital velocities closer to the Sun.
We calculated T coll (SE86) for the same meteoroids shown in Fig. 10B (a = 1.0 au, D = 50 µm) using the description given in Pokorný et al. (2018) and using the collisional lifetime multiplier F coll = 20 that was determined in Pokorný et al. (2014) for HTC meteoroids and subsequently used in numerous works.In Figure 10D, we show the ratio of our model and T coll (SE86) collisional lifetimes R = T coll /T coll (SE86).While there is some common ground (hashed areas in Fig. 10D show 0.56 < R < 1.78), our model predicts ∼10× longer collisional lifetimes for low-e and low-I meteoroids, while for meteoroids on high-e orbits our model predicts 1-2 orders of magnitude shorter collisional lifetimes than those estimated by T coll (SE86).These discrepancies are a natural consequence of assumptions made to calculate T coll (SE86).Steel & Elford (1986) assume the relative impact velocity is constant, which leads to an underestimation of destructive collisions near the Sun while overestimating the destructive collision rates for meteoroids on similar orbits.Grun et al. (1985) does not take into account the particle inclination and makes various assumptions for the velocity and density scaling in the solar system.In Figure 9E, we show that the Grun et al. (1985) and our model agree in predicted fluxes at 1 au for particles with D > 80 µm but show 1-2 orders of magnitude difference for smaller grains.Additionally, in Grun et al. (1985) the relative impact velocity of dust particles in the Zodiacal Cloud scales only with the heliocentric distance and ignores particle eccentricity and inclinations.This disagreement in predicted flux at 1 au is likely caused by a more complex size-frequency distribution at the source than the one assumed here.If we introduced a broken power law with a break point below D = 30 µm as suggested by Fixsen & Dwek (2002) and assumed a shallower power law index, this discrepancy would decrease.
Let us ignore for the moment the different flux and velocity distributions seen in our models and Grun et al. (1985).Eq. 4 indicates that Grun et al. (1985) assumes A s ≈ 300 J kg −1 , which corresponds to a value for our model parameter of P As = 3.Our JFC model (both scenarios) requires P As ≈ 5000.In other words, we require meteoroid binding energies ∼1000× larger than those assumed in Grun et al. (1985).
Figure 11 shows the model results assuming P As ≈ 10, (A s = 1000 J kg −1 ), a value comparable to that used in Grun et al. (1985).Since particles are much more prone to catastrophic destruction from much smaller particles than in our best fit model in both Scenarios, the larger meteoroids that are detectable by CMOR are efficiently removed before they can assume orbits with lower eccentricities and semimajor axis close to 1 au (Fig. 11A,C).The semimajor axis distribution of model meteors now peaks around a = 2.5 au, i.e. inside the main-belt, which leads to large discrepancy with the meteor radar observation and low goodness-of-fit values.The eccentricities of meteoroids are kept close to those of the source population and cannot reproduce the abundance of radar meteors with e < 0.8.This also impacts the observed meteor velocities, where the impacts with v ∞ < 30 km s −1 are missing from the model.Subsequently, the size-frequency distribution at 1 au is affected, mostly for meteoroids with diameters 100 < D < 1000 µm that are removed from the simulation before they can reach Earth.The larger meteoroids are not as affected because they do have shorter dynamical lifetimes due to their interaction with Jupiter (Yang & Ishiguro 2018).Frequent collisions also remove 63% of the mass flux at Earth compared to the best model fit from Scenario B. This result is not new and surprising.Nesvorný et al. (2011a) showed that using the original Grun et al. (1985) assumptions leads to early removal of meteoroids from the simulation and the fit to radar data cannot be achieved.Similar conclusions have been shown in Soja et al. (2019).These two works used a collisional lifetime multiplier, F coll , that increases the collisional lifetimes in the model by a factor of F coll .Nesvorný et al. (2011a) requires F coll > 30 for JFC meteoroids to provide a good match to their constraints, while Soja et al. (2019) requires F coll = 50 for meteoroids with D > 250 µm and F coll = 1 for smaller meteoroids.The reality is more complex as shown in Figure 12, where we show the values of F coll calculated from the ratio of the collisional lifetimes in our best model fit in Scenario B (A s = 500, 000 J kg −1 ) and approximate Grun et al. (1985) A s = 1000 J kg −1 for three different values of semimajor axis a and the full range of bound orbits in eccentricity and inclination assuming meteoroid diameter D = 50 µm.In Figure 12, we see that F coll varies by an order of magnitude for particles on various orbits, where F coll decreases with increasing eccentricity and inclination.In general, we can state that 30 < F coll < 100 for orbits with e < 0.8 for a wide range of semimajor axes in the inner solar system, which is in agreement with previous studies.However, seeing the complex picture resulting from our collisional grooming algorithm, we note that using simplistic estimates from analytically prescribed collisional lifetimes can under/over-estimate the collisional lifetimes by a factor of 3-5.7.1.JFC mass production rate and total mass loss rate and comparison to PSP One of the output of the Parker Solar Probe voyage to the outskirts of the Sun was an estimate of β-meteoroid production rates and conversely the total mass production rate in the inner solar system (Szalay et al. 2021).In this Section, we calculate the mass loss our JFC dust cloud suffers in collisions, compare the mass loss and its scaling with the heliocentric distance to Szalay et al. (2021), and estimate the JFC mass production rate needed to sustain the JFC portion of the Zodiacal Cloud.
To calculate the mass of dust that JFCs need to produce to sustain our best model fits in Scenarios A and B, we employ the following method.For each meteoroid diameter D and initial pericenter q ⋆ in our model, we calculate the total mass of meteoroids M nocoll (D, q ⋆ ) for a collision-less seed model and the collision-less dynamical lifetime T dyn (D, q ⋆ ).Using these two numbers we obtain the meteoroid production rate In Scenario A, we estimate m prod = 9, 360 kg s −1 while for Scenario B we obtain m prod = 5, 170 kg s −1 .Both of these values agree with Nesvorný et al. (2011a) who estimated 10 3 < m prod (JFC) < 10 4 kg s −1 .Our JFC production rates are approximately 5× smaller than those reported in Yang & Ishiguro (2018) who estimated 29, 000 < m prod (JFC) < 53, 000 kg s −1 .Rigley & Wyatt (2021) Figure 11.The same as Fig. 6 but now for Scenario B assuming meteoroid binding energy close to that used in Grun et al. (1985), i.e.P As = 10 (A s = 1000 J kg −1 ).Assuming much smaller A s than those providing best fit to our constraints leads to short dynamical lifetimes that prevent meteoroids detectable by CMOR to assume more circular orbits at semimajor axes close to 1 au.The size-frequency distribution at 1 au is also affected by removal of meteoroids with diameters 100 < D < 1000 µm.The mass accretion rate is at 37% of the best fit value in Fig. 6.
Figure 12.Estimated value of the collisional lifetime multiplier F coll for a particle with D = 50 µm and a full range of bound orbits phase space.We use Grun et al. (1985) binding energy estimate A s = 300 J kg −1 and A s = 500, 000 J kg −1 derived from our best fit model to derive the value of F coll .For the most of the orbital phase space the Grun et al. (1985) binding energy underestimates the collisional lifetimes by a factor of 40 − 100.
estimated the mass production rates of cometary fragmentation to the Zodiacal Cloud to have the mean value of m prod = 990 kg s −1 and the median value of m prod = 300 kg s −1 that are values an order of magnitude lower than our estimates.However, the same authors also find two epochs in their time dependent model where m prod = 6, 240 kg s −1 and m prod = 11, 110 kg s −1 , which is very close to our values.Additionally, using the Helios spacecraft data, Leinert et al. (1983) found m prod = 600 − 1000 kg s −1 , which is more similar to the values found in Rigley & Wyatt (2021) and lower than our estimates.
We can also calculate how much mass is produced in catastrophic collisions per unit of time by tracking the mass lost during each collision in our model.This is simply done by tracking the particle weights during the collisional grooming phase and multiplying by the mass that each particle in the model represents.In our model the, catastrophic collisions remove m loss = 430.6 kg s −1 and m loss = 164.5 kg s −1 for Scenarios A and B, respectively.This translates to 4.6% and 3.2% of the total mass of particles in the simulation.These estimated mass losses are close to the ∼100 kg s −1 estimated in Szalay et al. (2021) who used the β-meteoroid impact rates on the Parker Solar Probe to derive this number.The authors claim that 100 kg s −1 is a lower limit due to their assumption that all collisional products are β-meteoroids.This supports our findings since larger meteoroids (D > 1 mm) are likely to retain some mass in daughter products on bound orbits.Both of our scenarios allow a significant portion of mass lost in collisions to remain in bound orbits since our m loss is higher than the one reported in Szalay et al. (2021).The bound daughter collisional products would be then subject to PR drag and radiation pressure and experience similar fate as the parent dust grains of their sizes.
In contrast, Grun et al. (1985) estimate the m loss = 9000 kg s −1 for collisions inside 1 au.This value is similar to our meteoroid mass production rate that is needed to sustain the JFC meteoroid population in the Zodiacal Cloud.This would mean that all JFC particles in our model would be lost in collisions and no particles are lost to dynamical processes such as ejection of particles from the solar system, impacting one of the planets, and most importantly the particle disintegration near the Sun.Grun et al. (1985) estimates that 2.8% of particles are lost to dynamical processes (PR drag) and the rest is lost to collisions.Our model predicts that 95.4% and 96.8% of mass is lost to dynamical processes for Scenarios A and B, respectively.
If we decrease the particle strength by a factor of ∼1000 to A s = 300 J kg −1 which we used to estimate the Grun et al. (1985) collisional lifetimes, the meteoroid loss rate in collisions is m loss = 1647 kg s −1 and m loss = 787 kg s −1 for Scenarios A and B, respectively.This translates to 17.6% and 15.2% of particles lost due to collisions for Scenarios A and B, respectively.Even in such a case, the majority of mass is lost to dynamical processes since especially smaller particles are either blown out of the solar system or destroyed by spiraling too close to the Sun.The difference of our results to those reported in Grun et al. (1985) stems for a more comprehensive approach of our model by assuming the three-dimensional distribution of meteoroids in the Zodiacal Cloud and the relative impact velocities for each inter-particle collision.As we show in Figure 10, different combinations of orbital elements for particles with the same semimajor axis can lead to several orders of magnitude differences in collisional lifetimes.Also note that the A s = 300 J kg −1 model does not fit most of our observational constrains and our model requires values around A s = 500000 J kg −1 .
Our model also allows us to calculate the spatial distribution of mass produced in collisions per unit time, M + .Figure 13 shows the heliocentric distance scaling of M + for both scenarios A and  2021).The four realizations of our collisional model for JFCs are shown in color-coded solid lines: our preferred model A s = 5 × 10 5 J kg −1 (blue), 3× stronger particles A s = 1.5 × 10 6 (green) J kg −1 , 200× stronger particles A s = 10 8 J kg −1 (yellow), and 1000× weaker particles A s = 300 J kg −1 (grey).We also show the theoretically predicted scaling of M + based on brightness observations of the Zodiacal Cloud with a slope α = −3.33 (dashed black line).Our preferred model shows a good match to PSP measurements as well as the predicted slope down to heliocentric distance of 0.1 au, where the effects of the dust depletion zone near the Sun are evident.Panel B : The same as Panel A but now for Scenario B, where the total mass loss is m loss = 164.5 kg s −1 for the best fit parameters.In Scenario B, we assumed a smaller total cross-section of the Zodiacal Cloud and thus all models give smaller values of M + .Here, the 3× stronger particle model with A s = 1.5 × 10 6 J kg −1 underestimates M + compared to that derived from PSP measurements and provides an upper limit for A s in Scenario B.
The distribution of M + in the ecliptic plane is azimuthally symmetric as expected from our JFC simulation setup where we record particle and planet positions every 100 years and planets are sampled at random positions along their orbits (Fig. 14A).We can also calculate the scaling of M + with heliocentric distance, where get a value of M + ∝ r −3.83 hel between 0.1 < r hel < 1.0 au with a steeper decrease beyond 1 au M + ∝ r −5.5 hel .The slope of the M + scaling in the ecliptic is steeper (by 0.5 units) than the value for the volume average shown in Fig. 13.This is due to higher values of M + as well as the particle density close to the ecliptic, which is obvious from the side view shown in Fig. 14B.We fit the latitudinal profiles of M + using a function M + (β) = A exp(−B|z/r hel |) and obtained the value of the fitting parameter B = 8.3 ± 0.3, which is smaller than B = 2.1 derived by Leinert et al. (1981) from the Helios observations of the brightness of the inner Zodiacal Cloud.This shows that JFC meteoroids provide the close to the ecliptic component and the more diffuse component from Halley-type and Oort Cloud comets is needed to capture the entire image of the Zodiacal Cloud.

CONSEQUENCES FOR SOLAR AND EXO-SOLAR DUST AND DEBRIS CLOUDS
Our models show that the collisional strength of particles originating from Jupiter-family comets is ∼3 orders of magnitude higher than values commonly assumed in various solar system and debris disk studies.For dust originating from the Edgeworth-Kuiper Belt (EKB) objects, Kuchner & Stark (2010) assumed A s = 20 J kg −1 , whereas in Poppe (2016) and Poppe et al. (2019) the authors assumed that any collision results in a catastrophic disruption; i.e A s = 0.These values are > 25, 000× smaller than those we derived from our models and are even smaller than those assumed in Grun et al. (1985).It is understandable that EKB particles have a significant icy component that is likely to be less resistant to collisions and our JFC model cannot test that.The impact experiments of ice-silicate mixtures show that the higher percentage of ice content leads to higher impact yield and thus weaker particle binding energy (Koschny & Grün 2001).However, these results do not show 4 orders of magnitude changes in any measured quantity.Moreover, Benz & Asphaug (1999) show in their work that basalt and ice have similar collisional strength.For these reasons, we think that the effects of collisional grooming in Kuchner & Stark (2010), Poppe (2016), andPoppe (2019) are overestimated and should be revisited with a particle setup more resilient to collisions.As we show in Sec.7.1 the value of A s strongly influences the mass lost in collisions and the production of β-meteoroids, where more resilient particles can lead to 1-2 orders of magnitude smaller mass losses.This might lead to overestimation of the mass produced in collisions, shortening of the dynamical lifetime of particles in the model that are not able to get close to the host star, as well as underestimating the total mass of the dust source population if the total brightness of the dust cloud is used to constrain it.
Most exo-solar debris disk observations are modeled using three different codes that estimate the effects of collisions: ACE (Krivov et al. 2005), DyCoSS (Thebault et al. 2012), and LIDT-DD (Kral et al. 2013).We find that the values of A s or the critical specific energy Q * D used in models using these algorithms are mostly derived from the Benz & Asphaug (1999) work or are assumed to be smaller than in Benz & Asphaug (1999).In Benz & Asphaug (1999) the authors used smooth particle hydrodynamics (SPH) method to simulate collisions of cm-to-km sized bodies using basaltic and icy particles, where the impact velocities were considered to be between 0.5 km s −1 and 5.0 km s −1 .For basaltic material, the authors find A s = 9000 J kg −1 and A s = 3500 J kg −1 for impact speeds of 5 km s −1 and 3 km s −1 respectively.For icy material the values are: A s = 1600 J kg −1 and A s = 7000 J kg −1 for impact velocities 3 km s −1 and 0.5 km s −1 , respectively.In general, these values are ∼ 100× smaller than those derived in our work.Unfortunately, the Benz & Asphaug (1999) work does not discuss effects on particles smaller than 1 cm in diameter.This might be a potential source of discrepancy between our results and the values reported in Benz & Asphaug (1999) that are used in many debris disk studies.Such smaller values of A s lead, in general, to shorter dynamical lifetimes of all meteoroids in debris disk systems, higher mass production in collisions, and underestimation of total dust mass present in debris disks.However, each exo-solar system is different and the particular consequences of smaller/larger particle critical specific energy Q * D , that is linearly proportional to A s , need to be resolved by extensive modeling efforts that are beyond the scope of this article.

DISCUSSION
We provided the first implementation of a dynamical model for JFC meteoroids that includes iterative collisional grooming and which fits multiple observational constraints available in the inner solar system.Some aspects of the model require further improvement.We discuss these aspects here.
In this work, we assume that JFC particles dominate the inner solar system Zodiacal Cloud in terms of particle density, cross-section, and mass.Many studies support this picture (Nesvorný et al. 2010;Poppe 2016;Soja et al. 2019;Yang & Ishiguro 2018) and our work shows that JFC dust production can well reproduce all constraints we present in this article.By adding additional major dust populations such as main belt asteroid meteoroids, Halley-type comet meteoroids, and Oort cloud comet meteoroids, we do not expect significant changes in the best parameter fit we find for our JFC model.However, we expect that addition of particles with a broader inclination range, including some on retrograde orbits, will increase the rates of collisions in the innermost regions of the solar system because of how the particle density scales with the heliocentric distance (Poppe 2016;Pokorný & Kuchner 2019).
During our search for the best fit of our JFC meteoroid model to our set of constraints we explored two different scenarios with differing total particle cross-sections: (A) Σ ZC = 2.0 ± 0.5 × 10 17 m 2 based on Nesvorný et al. (2011a), and (B) Σ ZC = 1.12 × 10 17 m 2 based on Gaidos (1999).We showed that we were able to reconcile the Scenario A fit by using the mass accretion rate of JFCs at Earth from Carrillo-Sánchez et al. (2016), whereas in Scenario B we show that the Gaidos (1999) estimate of the total Zodiacal Cloud cross-section matches well the mass accretion rates from Carrillo-Sánchez et al. (2020).We conclude that more constraints are needed to decide which of the two presented scenarios truly represents the current state of the Zodiacal Cloud.
Our collisional grooming model itself is still an approximation of reality where we assume that meteoroid collisions have either no effect or they catastrophically disrupt particles involved in the collision.We also assume that the dust produced in collisions does not further influence the collisional evolution of dust in our model and that all collisional fragments are beta-meteoroids.We also assume that the Zodiacal Cloud and our Jupiter-family comet meteoroid model is unchanging in time and is averaged over timescales longer than orbital periods of the eight planets.This is not entirely correct since observations and models show that the Zodiacal Cloud experiences time variability on thousand to million year timescales (Napier 2001;Koschny et al. 2019;Rigley & Wyatt 2021).
While we attempted to acquire a very broad selection of model constraints, we acknowledge that more data sets are available that were not discussed in detail in this article.We are seeking to adopt an even broader set of observations and measurements that would allow us to better constrain the meteoroid production rates in collisions, provide independent measurements for mass accretion rates at various planets, and set constraints on the global shape of the Zodiacal Cloud for a broad range of dust and meteoroid sizes.With our model available to the general public, we aim to keep the model updated via community updates on the project website.
One of the potential drawbacks of our current model for JFC dust population is the assumption of a single power law of particles generated by the source population.We decided to avoid a more complicated SFDs such as the broken power law due to the additional fitting parameters (2 for the broken power law).This increase in the size of the parameter space would significantly extend the total processing time required to find the best fits to our suite of constraints.Introducing a broken power law would influence the particle number fluxes shown in Figs. 6 and 9 as well as the particle density with respect to heliocentric distance.It would also change the collisional lifetimes of particles with diameters close to and below the break point of the power law.We do not expect significant changes for the strength of the particles A s derived in this work since the most collisionally impacted are larger particles with D > 100 µm and we expect the break points to be around D ∼ 30 µm based on our initial analysis and the FIRAS measurements of the Zodiacal Cloud spectrum by Fixsen & Dwek (2002).

Time and memory consumption of the collisional grooming process
The collisional grooming process run time is influenced by two major factors: (1) the total number of records in the model/cloud and (2) the number of iterations needed to attain a steady state.Most of our modeling runs took 5-6 iterations to reach the first steady state and then 1-2 further iterations to reach the second, rescaled steady state.The iteration time varied between different computer setups, ranging between 3000 and 4000 seconds per iteration, or approximately 6-8 hours for one model run.
The original memory required per run was over 128GB RAM which limits the model deployment on most computers available to us.Therefore, we implemented several tweaks that reduce the memory requirements without decreasing the precision of the collisional grooming method.This was done through the aggregation of particles inside each volume element that had the same diameter and similar velocity vectors.The particle position information is not retained once the particle is assigned to a particular volume element, thus it plays no role in particle aggregation.The velocity vector similarity is checked in two steps: first the angular difference ⃗ v 1 • ⃗ v 2 /(∥v 1 ∥∥v 2 ∥) = cos ϕ is checked, and then the velocity vector size ratio is checked R v = | log ∥v 1 ∥ − log ∥v 2 ∥|.In our model, we aggregate particles with the same D that have a velocity vector angular difference < 10 • (cos ϕ > 0.985) and their vector size ratio is < 10% (R v < 0.0414).Making this approximation reduced the memory usage to approximately 17.3 GB of RAM.This memory requirement still limits deployment on base-model laptops, which usually have 8-16 GB of RAM, but it suits the capabilities of a typical "professional" desktop computer with 32+ GB of RAM.We tested the loss of precision

Figure 1 .
Figure 1.Maximum target and impactor size ratio for destructive collisions R max .If the target is larger than R max then the target is not destroyed during the collision.This Figure shows values of R max for four different values of A s .In Panel A we show A s = 10 2 J kg −1 , which is equal to A s = 3.02 × 10 6 erg g −1 used in Krivov et al. (2005).Panels B-D show R max for A s = 3 × 10 4 , 3 × 10 5 , and 10 6 J kg −1 which are cases of interest in this manuscript.

Figure 3 .
Figure 3. Initial conditions for JFCs in our numerical models.Panel A shows the cumulative distributions of semimajor axis a for 10 different values of critical perihelion distance q ⋆ (color-coded solid and dashed lines) and the currently observed distribution of JFCs with comet total magnitude parameter M 1 ≤ 12. Panels B-D are the same as panel A but show instead the inclination (panel B), eccentricity (panel C), and perihelion distance (panel D) cumulative distributions.The legend also shows the total number of bodies for each perihelion distance set, N .The data for currently observed JFCs were obtained on May 20th, 2021 from https://ssd.jpl.nasa.gov/sbdbquery.cgi.ASCII version of the cumulative plots is available at the project GitHub.

Figure 4 .
Figure 4. Goodness-of-fit function variations for different values of model parameters P. Panel (A) shows the influence of P α on the goodness-of-fit functions for fixed P γ = 0.15, P As = 5500, P Bs = −0.24.Goodness-of-fit functions are color-coded: χ 2 R is shown in yellow, χ 2 Ṁ in pink, χ 2 O in blue, and χ 2 α in black.The gray region shows the parameter space we used for our focused fitting P α ∈ [4.2, 4.6].Panels (B-D) show the same as panel (A) but for parameters P As , P γ , P Bs , respectively.The focused fitting parameter space is (B) P As ∈ [1000, 8000], (C) P γ ∈ [0.125, 0.175], and (D) P Bs = −0.24.

Figure 5 .
Figure 5. Variations of the goodness-of-fit function χ 2 tot of our Scenario A JFC model for 216 different parameters value combinations (P α ,P As ,P γ ).The three symbols represent three different values of P γ : (×; 0.125), (•; 0.15), (+; 0.175).For clarity we offset × and + symbols from • but all three symbols have the same values of (P α ,P As ).Each of the model realizations is color-coded by the value χ 2 tot .The region of the best fitting parameters is denoted by the black rectangle: 4.325 ≤ P α ≤ 4.425 and 1500 < P As < 8500.Panel (A) shows the full range of the total χ 2 tot where we recorded the minimum value χ 2 tot = 1.13 for (P α ,P As ,P Bs , P γ ) = (4.40,8000, −0.24, 0.125).Panel (B) contains the same information, but the color-scale now spans those values of χ 2 tot within 10% of the minimum value (1.13).Panel (C) shows the variation in the goodness-of-fit for mass accretion at Earth, χ 2 Ṁ , while panel (D) shows the variation in the goodness-of-fit for size-frequency distribution slope, χ 2 α .

Figure 6 .
Figure 6.Overview of the best model for scenario A and its ability to fit our constraints.Panels (A-D) show the semimajor axis (A), inclination (B), eccentricity (C), and hyperbolic excess velocity (D) as observed by CMOR (orange solid line) and as reproduced in our model (black solid line).The gray solid lines show the collisional grooming iterations, where the lighter gray tones show the early iterations and the darker show more final ones.The blue line represents the model residuals.Panel (E) shows the size-frequency distribution of particles impacting Earth (black solid line), the broken-power law fit (orange dashed line) and Grun et al. (1985) interplanetary (pink dashed line) and lunar (blue dashed line) particle fluxes.The Grun et al. (1985) fluxes are multiplied by a factor of 3 to account for gravitational focusing effect that Earth experiences with respect to the Moon and interplanetary space (Pokorný et al. 2018).The legend shows broken power law fit parameters and the fitting function.All fluxes are cumulative.Panel (F) shows the distribution of particle density with heliocentric distance measured in the ecliptic plane for our model (black solid line) and the fitted exponential (orange dashed line).The fit parameters and fitting function are shown in the legend.Finally, the blue panel shows the value of M Ṁ and the goodness-of-fit function values.The red panel shows the model parameters P.

Figure 8 .
Figure 8.The same as Fig. 5 but now for Scenario B using the Gaidos (1999) total Zodiacal Cloud cross-section Σ ZC = 1.12 × 10 17 m 2 for our JFC model.The three symbols represent three different values of P γ : (×; 0.175), (•; 0.20), (+; 0.225).For clarity, we offset × and + symbols from • but all three symbols have the same values of (P α ,P As ).Each of the model realizations is color-coded by the value χ 2 tot .The region of the best fitting parameters is denoted by the black rectangle: 4.1725 ≤ P α ≤ 4.225 and 1500 < P As < 8500.Panel (A) shows the full range of the total χ 2 tot where we recorded the minimum value χ 2 tot = 0.299 for (P α ,P As ,P Bs , P γ ) = (4.20,5000, −0.24, 0.20).Panel (B) shows the values within 10% of the minimum value χ 2 tot = 0.292.Panel (C) shows the variations of the goodness-of-fit for mass accretion at Earth χ 2 Ṁ , while panel (D) shows the same for the goodness-of-fit χ 2 α for size-frequency distribution slope α.

Figure 9 .
Figure 9.The same as Fig. 6 but now for the best model for Scenario B.

Figure 10 .
Figure 10.Scenario B collisional lifetimes for D = 50 µm meteoroid.Panels (A)-(C) show the collisional lifetime in years for the full range of bound orbits; eccentricity 0 < e < 1, inclination 0 < I < 180 • where the color coding represents logarithmically spaced collisional lifetime.Panel (A) is for a particle with semimajor axis a = 0.3 au, panel (B) for a = 1.0 au, and panel (C) for a = 3.0 au.Contours show boundaries between color-coded logarithmic steps.Panel (D)  shows the ratio between collisional lifetimes in our model T coll and T coll (SE86) with F coll = 20 used in e.g.,Pokorný et al. (2014Pokorný et al. ( , 2018)), i.e. the color-coding shows R = T coll /T coll (SE86).The hatching shows areas within a factor of 1.78 of unity, i.e R = 0.56 − 1.78.Our model provides ∼3× longer collisional lifetimes for meteoroids on low-e and low-I orbits and 1-2 orders of magnitude shorter collisional lifetimes for high-e and high-I orbits (dark blue regions).

Figure 13 .
Figure 13.Panel A:The heliocentric distance scaling of the collisional mass production rate M + for Scenario A, where the total m loss = 430.6 kg s −1 for the best fit parameters.The dashed orange line shows the PSP measurements derived fromSzalay et al. (2021).The four realizations of our collisional model for JFCs are shown in color-coded solid lines: our preferred model A s = 5 × 10 5 J kg −1 (blue), 3× stronger particles A s = 1.5 × 10 6 (green) J kg −1 , 200× stronger particles A s = 10 8 J kg −1 (yellow), and 1000× weaker particles A s = 300 J kg −1 (grey).We also show the theoretically predicted scaling of M + based on brightness observations of the Zodiacal Cloud with a slope α = −3.33 (dashed black line).Our preferred model shows a good match to PSP measurements as well as the predicted slope down to heliocentric distance of 0.1 au, where the effects of the dust depletion zone near the Sun are evident.Panel B : The same as Panel A but now for Scenario B, where the total mass loss is m loss = 164.5 kg s −1 for the best fit parameters.In Scenario B, we assumed a smaller total cross-section of the Zodiacal Cloud and thus all models give smaller values of M + .Here, the 3× stronger particle model with A s = 1.5 × 10 6 J kg −1 underestimates M + compared to that derived from PSP measurements and provides an upper limit for A s in Scenario B.

Figure 14 .
Figure 14.Panel A: Distribution of M + in the ecliptic plane in our preferred model in Scenario B. The plot shows the average M + within 0.025 au from the ecliptic plane.The value of M + in the ecliptic is radially symmetric in our model and decreases as M + ∝ r −3.83 hel between 0.1 < r hel < 1.0 au, which is steeper than the radially averaged value of M + shown in Fig. 13.The color-coded contours show the order of magnitude increase in M + where our model predicts 9 orders of magnitude difference between M + at r hel = 5 au and r hel = 0.1 au.Panel B: The same as Panel A but now for a slice through the Y -axis, a side view of the JFC dust cloud where the values of M + are averaged over 0.025 au away from the Y -axis in our model.Due to rotational symmetry these two images can be used to reconstruct the entire 3D image of the distribution of M + in our JFC model.

Table 1 .
Initial number of particles in each bin of our dynamical model of JFC meteoroids.