{'fraction_non_alphanumeric': 0.03479965432202854, 'fraction_numerical': 0.013932508483864849, 'mean_word_length': 4.212480773456384, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Two-body model for the spatial distribution of dust ejected from an atmosphereless body

We present a model for the configuration of noninteracting material that is ejected in a continuous manner from an atmosphereless gravitating body for a given distribution of sources. The model is applicable to material on bound or unbound trajectories and to steady and nonsteady modes of ejection. For a jet that is inclined to the surface normal, we related the distributions of ejection direction, velocity, and size to the phase-space number density at the distance from the source body. Integrating over velocity space, we obtained an expression from which we inferred the density, flux, or optical depth of the ejected material. As examples for the application of the code, we calculate profiles of the dust density in the Enceladus plume, the pattern of mass deposition rates around a plume on Europa, and images of optical depth following the nonstationary emission of material from a volcano on Io. We make the source code of a Fortran-95 implementation of the model freely available.


Introduction
The ejection of material from the surfaces of atmosphereless bodies is a ubiquitous phenomenon in the Solar System. Prominent examples are comets, active asteroids, ejecta clouds from hypervelocity impacts, or plumes from active satellites. For the dynamics of the ejected material, it is in many cases possible to neglect any other forces than the mass point gravity of the source body to a good degree of approximation. This is for instance the case for impact-generated dust clouds around planetary satellites as were detected around the Galilean moons (Krüger et al., 1999(Krüger et al., , 2003 or the Moon (Horányi et al., 2015), or dust plumes ejected from cryovolcanically active satellites (Spahn et al., 2006;Porco et al., 2006;Southworth et al., 2015). We note that for higher-order gravity terms to be negligible, the source body does not necessarily need to be spherical. For instance, mass-point gravity can be a good approximation to describe dust ejection from a satellite with surface topography.
In this paper we derive a semianalytical model to assess the spatial configuration of the emitted dust. The model relates the distribution of dust sources on the surface of the atmosphereless body and the parameters of ejection (e.g., source strength or directional and velocity distribution) to observable quantities such as number density, fluxes, or optical depth. The mathematical foundations are described in Sect. 2. Expanding on work in the literature Sremčević et al., 2003), our model can handle emission through inclined jets, and we develop a method for carrying out two of three integrations over the velocity distribution analytically. The code that implements the new model, carrying out the one remaining integration numerically, is called DUDI (for "dust distribution") and is freely available under the GNU General Public License on https://github.com/Veyza/dudi. Aspects of the numerical algorithm for the integration are outlined in Sect. 3. Examples for an application of the model to current problems in planetary science are given in Sect. 4, including cases of steady and nonsteady dust emission.
2 Mathematical formulation 2.1 Phase-space density We followed the derivations by Krivov et al. (2003), Sremčević et al. (2003), and Postberg et al. (2011) to relate the phase-space density of dust in a certain point of interest in space to the distributions that describe the ejection of the dust from a source on the moon surface. We generalized the existing model to allow emission from a point source in a direction that is not normal to the surface, with an axisymmetric distribution of ejection angles around this direction. Moreover, we allowed for a general coupling of the distribution of ejection velocities and grain size.
Our model was developed initially to fit in situ measurements by the Cassini Cosmic Dust Analyzer at the Saturn satellite Enceladus. For convenience, we use the words "spacecraft position" or "spacecraft coordinates" from now on to denote the point in space at which the dust density is calculated. We also use the term "density", which at any point can be understood as the number density, mass density, the average radius of dust particles, or the cross section that is covered by the dust at the spacecraft position. The model allows us to obtain any of these quantities with a change of only one parameter.
We first consider a stationary process. We can equate the differential number of dust particles in a certain point of phase space dn = n(r, α, β, v, θ, λ, R)r 2 sin αdrdαdβ v 2 sin θdvdθdλdR (1) to the number of particles ejected from the satellite surface dn = γdt f (α M , β M , u, ψ, λ M , R) sin α M dα M dβ M du sin ψdψdλ M dR.
For the two-body problem, the Jacobian can be obtained analytically (see Sremčević et al. (2003)), We assume that the distribution function f factorizes, and that the distributions of the source position (f α M ,β M ), ejection speed (f u ), ejection direction (f ψ,λ M ), and size of the ejected dust particles (f R ) can be defined separately, It is physically plausible that the distribution of the ejection speed of the dust particles depends on the grain size (Schmidt et al., 2008;Postberg et al., 2011), which we emphasize with the notation f u (u, R).
We describe the position of the point source located at the coordinates (α 0 M , β 0 M ) on the surface of the spherical moon as the product of two Dirac δ-functions, The term sin α M in the denominator comes from the normalization. To formulate the directional distribution of ψ and λ M so that it describes the ejection of dust that is axisymmetric around the axis of an inclined jet, we consider two coordinate systems centered at the location of the point source. The Z -axis of system (X, Y, Z) points along the local normal to the surface, and the X-axis points toward the local north. TheZ-axis of system (X,Ỹ ,Z) is aligned with the axis of the jet. The axisX lies on the line of nodes, so that the angle η * measured from X toX is related to the jet azimuth η as η * = η − π/2. We define azimuth angles always clockwise from the local north, allowing a direct comparison to the derivations in Krivov et al. (2003) and Sremčević et al. (2003). Then, the transformation of the (X, Y, Z) coordinate system Particle speed at the spacecraft position u Particle speed at the moment of ejection θ Angle between the particle velocity and position vectors ψ Initial angle between the particle velocity and position vectors λ Azimuth angle of the particle velocity, measured clockwise from local north.
λ M Azimuth angle of the particle initial velocity, measured clockwise from local north. R Radius of the particle γ Rate of dust particle production n Phase-space density of particles with fixed radius f Distribution describing the dust ejection process ζ Zenith angle of the source symmetry axis η Azimuth of the source symmetry axis measured clockwise from local north η * Auxiliary angle used in the derivation of the expression for the ejection direction distribution in case of a tilted symmetry axis to the (X,Ỹ ,Z) coordinate system may be performed as two subsequent rotations, as shown in Fig. 2. The first rotation is clockwise around the Z-axis with angle η * . The second rotation is counterclockwise around theX-axis with angle ζ. Let (ψ, λ M ) and (ψ,λ M ) be the polar angle and azimuth in the two systems (X, Y, Z) and (X,Ỹ ,Z), respectively. The distribution of the ejection direction we wish to use can be formulated in a simple manner in terms of the variables (ψ,λ M ) because the distribution is symmetrical with respect to the axisZ. However, in Eq. (3) and in the formulae that are to be derived later in the course of solving the two-body problem, we have to deal with the angles (ψ, λ M ) that are defined in the local horizontal coordinate system. Therefore a replacement of the coordinates must be performed in the expression to obtain the distribution of ψ and λ M , which can be used in further calculations. The desired function f ψ,λ M (ψ, λ M ) can be obtained by multiplication with the Jacobian, To expressψ andλ M through ψ and λ M , we consider a unit vector k pointing in an arbitrary direction (Fig. 2). In both coordinate systems the vector can be defined by Cartesian coordinates related to the corresponding polar coordinates as (we recall that λ M andλ M are azimuthal angles counted clockwise from the X andX axes, respectively) Figure 2: Two coordinate systems centered at the location of the dust source on the surface of a spherical body. The Z-axis of system (X, Y, Z) is normal to the surface, and X points to the local north. TheZ-axis of the coordinate system (X,Ỹ ,Z) is aligned with a jet that is tilted from the surface normal by an angle ζ.
The Cartesian coordinates of k in the two systems are related through the rotation matrix, which can be expressed in terms of η, Using Eqs. 8 and 9, we obtaiñ ψ = arccos(cos ζ cos ψ + cos(η − λ M ) sin ζ sin ψ), This gives the Jacobian

Integration
To compute the density of dust at the point (r, α, β), we must integrate Eq. (3) over all possible velocities and over all possible particle sizes, Here, is defined in a similar way as in Postberg et al. (2011). The parameter p defines the moment of the size distribution related to the quantity we are interested in. Using p=0, we obtain the number density of particles in the specified range of sizes. Setting p=1 gives the average radius of the grains per unit volume. For p=2 we obtain the average cross section of the dust particles per volume. This setting is used below to compute the geometrical optical depth of the dust population. Finally, p=3 gives the average volume occupied by dust grains per unit volume. This setting is used to compute the mass density of the dust. For more details of the evaluation of G p u (R min , R max ), see Appendix B. We replace variables in the argument of the δ function in equation (13) as to integrate over θ and λ analytically. Eq. (15) is derived in greater detail in Appendix A. Here, F (θ, λ) represents the integrand of equation 13, while θ i and λ i are the roots of the equation All the necessary dependencies between the variables in question can be found from spherical trigonometry Sremčević et al., 2003), for instance, α M = arccos (cos α cos ∆φ(θ) − sin α sin ∆φ(θ) cos λ) and The spherical triangle used to obtain these relations is shown in Fig. 1. The angle ∆φ is the angle between the position vectors of the spacecraft and the source location on the moon. Because θ enters expressions 17 and 18 only through ∆φ, the partial derivatives of α M and β M with respect to θ can be computed as partial derivatives with respect to ∆φ multiplied by ∂∆φ/∂θ. The Jacobian reads and our final formula is The integration over v must be carried out numerically. The lower integration limit is restricted by the minimum energy, or minimum semimajor axis, of the orbits that pass through the two points (r M , α M , β M ) and (r, α, β). It is obtained from where a min = r + r M 4 + 1 2 The upper limit is constrained by the maximum ejection speed, The angle ∆φ in the expressions above is the angle between the position vector of the dust source, r M , and the position vector of the spacecraft, r. For a fixed position of the source and spacecraft, the angle ∆φ is also fixed, but it formally depends on v and θ. In the integrand of equation 20, the value v is given and ∆φ(θ) is a function of the variable θ alone. From the conservation equations of the two-body problem, we obtainp = 2r 2ṽ2 sin 2 θ, which give the relation between ∆φ and θ. Following the notational convention of Krivov et al. (2003), we use dimensionless variablesr = r/r M andṽ = v/v escape , where v escape is the escape velocity on the satellite surface. The angles φ and φ M are the true anomalies at r and r M , respectively. Equations (24) -(28) can be used to calculate the partial derivative ∂∆φ/∂θ. However, it is not possible to reverse these expressions to obtain θ from a given value of ∆φ analytically. The desired θ * i are the values of θ that (for a given v) lead to a ∆φ satisfying equations 17 and 18. We use a geometrical approach to calculate all possible θ * i . We consider the two-body problem, therefore the motion is restricted to a plane. We define a two-dimensional coordinate system in the plane containing the vectors r and r M . The origin is located at the moon center and r points along the x-axis. We know the lengths of r and r M as well as the angle ∆φ between them. Then the coordinates of the vectors r and r M in the plane are (r, 0) and (r M cos ∆φ, r M sin ∆φ).
We wish to obtain θ , which is the angle between r and the grain velocity vector, which is tangential to the particle trajectory at the point r. This angle can be calculated when we know the equation of the trajectory, which is either an ellipse or a hyperbola. At this step (calculating the integrand of equation 20), we know the value of the particle speed and the distance to the moon center, which determines the orbital energy, and thus, the semimajor axis a. We also know whether the particle moves along an ellipse (negative orbital engery) or along a hyperbola (positive orbital energy).
We first consider the case of an ellipse (Fig. 3). One of the focal points (F 1 ) is located in the origin, which is the center of the moon. To draw the ellipse, we must find the position of the second focus. At each point of the ellipse, the sum of the distances to the focal points is a constant equal to 2a. Therefore the second focus must be removed from point r by 2a − r and from point r M by 2a − r M . These conditions are met at the intersection points of two circles centered at r and r M with radii 2a − r and 2a − r M , respectively. This leads to the conditions By solving the system of equations 29, we find two possible positions for the second focal point, F 1 2 and F 2 2 . In either case, we can calculate the eccentricity e of the ellipse, the true anomaly at point r, and, finally, the two solutions for θ using Eqs. (30) -(33). Knowing cos φ M is sufficient to obtain φ M because we know that the ejection point r M cannot have a true anomaly φ M > π. We have These equations determine the solutions θ * i used in equation 20 if the particle travels from r M to r along the shorter arc of the ellipse. However, there are cases when a particle reaches r over an arc of 2π − ∆φ (Fig. 5), leaving r M in the opposite direction. To distinguish this case, we must recalculate r from the obtained value for φ using and verify that it matches the starting value for r that we used to obtain φ. If this is not the case, then ∆φ in Eq. (32) must be replaced by 2π − ∆φ, which corresponds to the motion along the same ellipse, but in the opposite direction. This case is relatively rare, and in the examples we explored was only encountered at large distances from the source. In case of hyperbolic motion (a < 0, Fig. 4), we follow almost the same line. For every point on a hyperbola, the difference between the distances to the focal points is the same. Because r M and r must be distances to the nearest focus, the system of equations for the coordinates of the second focal point reads Equations 30 -33 remain the same for the hyperbolic case. However, when the coordinates of the second focus are found, we must make sure that the particle does not pass the pericenter on its way from r M to r. For this purpose, we verify that the points r M and r lie on the same side of the line F 1 F 2 . If this condition is not satisfied, the solution is rejected. In Fig. 4 the hyperbola with second focus at the point F 2 2 does not meet this condition. Therefore only one hyperbolic trajectory is possible to get from r M to r. Furthermore, the motion along the hyperbola is possible only in one direction. The value for sin φ is always positive and ∆φ = φ − φ M .
The values of λ i can be inferred from spherical trigonometry. The two solutions can be either identical or they differ by 180 • because the motion is restricted to a plane, The sign in Eq. (37) depends on the specific orientation of r and r M relative to the direction of zero-longitude. When the particle travels from r M to r over the angle of 2π − ∆φ , the signs of sin λ and cos λ both change because the sign of the sin ∆φ term in the denominator changes.
As soon as the values of (v, θ * i , λ i ) that satisfy the orbital geometry are known, they may be used to calculate the corresponding (u, ψ i , λ M i ) from Eqs. (38) -(41) and the integrand in Eq. (20) is fully determined,

Nonstationary case
Nonstationary dust ejection can be modeled by allowing a time-dependent production rate γ in Eq. (20) that will result in a time-dependent spatial distribution of the dust. When we determine the orbital geometry for a fixed velocity value at the given point in space (Sect. 2.2), the time ∆t required for traveling from r M to r along the Keplerian orbit can be calculated from Kepler's equation. Thus, we know that the properties of the dust configuration at location r and time t are caused by the production of dust at the source location r M with the rate γ(t − ∆t). In this case, the production rate cannot be put outside the integral and Eq. (20) becomes The two solutions for ∆t i follow from Kepler's equation using the two solutions for eccentricity from equation 30 and computing the eccentric anomaly with the half-angle formula from the true anomaly given by equation 32.

Singularities of the coordinates
There are coordinates for which solutions of equations (20) or (42) cannot be obtained (see Table  2). The problem arises from the use of a spherical coordinate system. In practice, it is possible to avoid the singularities by carefully choosing the pole axis of the coordinate system after the source location and the detector position of interest are known. If coordinates close to the singularities need to be evaluated, then the stability and accuracy of the numerical integration becomes a challenge. However, double-precision calculations allow approaching the singular angles as close as Point source of dust is located at the north or south pole, α M = 180 • its latitude is not defined, and no spherical triangle from Fig. 1 exists α = 0 • Spacecraft is located directly above the north or south pole, α = 180 • its latitude is not defined, no spherical triangle from Fig. 1 A particle is ejected exactly along z = 0 • the jet axis, the azimuth of ejection is undefined, and the Jacobian in Eq. (12) diverges

Numerical integration
In this section we present and discuss the algorithm for the numerical solution of equations 20 and 42.
We have implemented this algorithm in a code written in Fortran-95, which we call dust distribution (DUDI). The source code with technical documentation and instructions for usage and compilation is freely available under the GNU General Public License on https://github.com/Veyza/dudi. The makefile provided for compilation uses the gfortran compiler. DUDI can be compiled and run without installing additional libraries. The library of OpenMP, which is used to speed the computation up, is included in the compiler. DUDI allows us to compute the number density of dust or related quantities at given points in space as a mean radius, average cross section, or mass density of the dust grains ejected from the surface of a spherical body without an atmosphere. The input data are the spacecraft coordinates, the properties of the source (i.e., the location and the distributions of the direction and speed of the ejection) and the three parameters of G p u . The first preliminary step is to calculate G p u on a grid of u-values. The array of pairs (u, G p u ) is later used to interpolate G p u (R min , R max ) for the actually required value u under the integrand of equations 20 or 42. The second preliminary step is to compute the values of ∆φ and ∆β for the given positions of the source and the spacecraft.
Then we proceed directly to the numerical evaluation of the integral over velocity in Eq (20) in the stationary case, or Eq. (42) if there is a time-dependence. The ejection speed distribution implies a certain lower limit for the possible ejection speed u min . At any point r, this minimum ejection velocity restricts the corresponding minimum velocity at spacecraft position v 0 min , which may be higher than the lower integration limit given by Eqs. (22) -(21). The actual numerical integration is performed over the interval where (v min , v max ) and (v 0 min , v max ) overlap. The minimum and maximum ejection speed u min and u max (in m/s) must be explicitly specified as a property of the dust source, along with the expression for the ejection speed distribution.
At each step of the integration, we find for a given v, r, and ∆φ the solutions for the angles θ and λ as described in Sect. 2.2. We control the accuracy of the solution for θ by recalculating the value of ∆φ from Eqs. (24)-(28) and comparing it to the starting value of ∆φ as given by the positions of the dust source and the spacecraft. We require the difference between the two values of ∆φ to be smaller than 10 −4 ∆φ. In most cases the accuracy is much better, but it may degrade for certain values of v near the singularities of the coordinates (see Sect. 2.4) or when r ≈ r M . Even for poor accuracy for θ , this means that the accuracy decreases only for one or two integration steps. The error in the final result is smaller. We consider the relative accuracy of 10 −4 sufficient for the θ solutions. If it is worse, a warning is produced by the program.
We divide the integration domain into two regions as (v min , v par ) and (v par , v max ) to obtain the number density of the particles on elliptic and hyperbolic (escaping) trajectories separately. Here v par (the subscript stands for "parabolic") is the minimum escape velocity at radial distance r.
When v min > v 0 min , the integrand in Eq. (20) has a pole at v = v min . We replace v min from Eq. (21) with v min + ∆, where ∆ = 10 −10 has turned out to be a good choice to evaluate the integrand near the pole to reasonable accuracy in a stable manner. To better resolve the pole, we additionally subdivide the elliptic part of the integral into two parts that are treated separately. The first integration subinterval is (v min , v 1 ), where v 1 = v min + 10 −4 (v par − v min ) is the domain that contains the pole. Here integration is performed using the trapezoidal rule with a large number of supports that become denser toward v min as where N is the number of supports. Then we use the Gauss-Legendre quadrature of a moderate order to compute the integral from v 1 to v par . The integration over the hyperbolic velocities is also performed with a Gauss-Legendre quadrature. The choice of the quadrature order depends on the integration domain and ejection speed distribution. The nodes and weights of the Gauss-Legendre formula are tabulated in our code for the following orders: 5, 10, 20, 30, 40, and 50. Figures 6 -9 show examples for the behavior of the integrand in the three domains. Depending on the choice of the ejection speed distribution f u (u, R), the integrand may decrease (Fig. 7) or increase (Fig. 8) toward higher velocities. Remarkably, the integrand can jump at v = v par (Fig.  9) if a significant part of the dust number density is due to the particles on their way back to the moon after passage of their apocenter, The sharpness of the pole varies. At each integration step, the given velocity v determines the values of θ that enter the integrand of equation 20 through the derivative ∂∆φ/∂θ and ψ(θ). The latter is needed to compute the ejection angle distribution f ψ,λ M (ψ, λ M ). The factor ∂∆φ/∂θ is the reason for the pole. Its value depends on θ and on the spacecraft position relative to the source. The pole is less strongly peaked if the value of ψ(θ) corresponds to a very unlikely ejection direction. Thus, the sharpness of the pole depends on the spacecraft position relative to the source position and also on the ejection angle distribution, and so does the number of integration steps required to achieve a given accuracy goal.
The accuracy can be estimated from Eq. (44), where P is the value of the pole integral (between v min and v 1 ) obtained with N steps for the integration with the trapezoidal rule, and N max is the maximum reasonable number of steps. N max is limited by accumulated rounding errors, and we determine its value in test integrations. I(N max ) is the sum of P (N max ) and the remaining part of the integral between v 1 and v max . In this way, we can quantify the discrepancy induced by the pole integration in the final result, We require ≤ 10 −3 and perform tests to determine the corresponding number of pole integration steps N necessary to achieve this goal. This number we compute for different spacecraft positions relative to the source and to the axis of ejection symmetry (the polar angle in the coordinate systemXỸZ in Fig. 2, in the following denoted by ξ). We adopt for the ejection direction distribution. Normalization in this expression does not matter for an evaluation of from Eq. (44). We vary the parameters ψ max and ω, along with the polar angle of the ejection symmetry axis, to investigate the behavior of the pole for different ejection distributions, of which two main classes can be defined. The "jets" are the sources with a preferred direction of ejection, and the "diffuse sources" have no such direction. We find that for diffuse ejection (ω = 45 • and ψ max = 45 • ), N = 15 is a sufficient number of supports to achieve ≤ 10 −3 at any spacecraft position. For a vertical jet (ψ max = 0 • and ω in the range of 3 • and 5 • ), the value of P can be neglected for all ξ > 40 • and N = 15 is sufficient for ξ < 40 • . However, for a narrow and inclined jet, we find that there are points where a large number of supports is required to integrate the pole accurately. The narrower and the more inclined the jet, the greater the number of these points and the greater the required N . We focus on the worst-case scenarios generally to constrain an optimal number of steps required for the pole integration.
Empirically, we find that an accuracy of 10 −3 can be achieved with a minimum number of steps when an exponent of k = 4 is used in Eq. (43). Figs. 10 and 11 show examples of how the minimum number of steps required to integrate the pole with the given accuracy is distributed over r and ξ values. N = 0 means that the pole does not have to be integrated at all because its value is negligible. For jets, we select the number of supports N for the pole integration based on the distributions shown in Figs. 10 and 11 as follows:

Applications
In this section we present three applications of the model to phenomena of scientific interest in the Solar System. The purpose of these examples is to demonstrate the wide range of applicability of the model. We leave a rigorous scientific analysis of these problems with a comprehensive comparison to data for future work.

Density profile of the Enceladus dust plume
On July 14, 2005, the Cassini spacecraft performed a flyby at the Saturnian moon Enceladus (labeled E2). During the flyby, the High Rate Detector (HRD), a subsystem of the Cassini Cosmic Dust Analyzer instrument (Srama et al., 2004), measured the number density of dust particles in the vicinity of the satellite. The significant increase in dust density near Enceladus (see Fig. 12), about one minute prior to the closest approach of the spacecraft to the satellite, was the first in situ measurement of particles in the Enceladus dust plume (Spahn et al., 2006). Dust and vapor are emitted from four fissures called the tiger stripes in the anomalously warm south polar terrain of Enceladus (Spencer et al., 2006;Porco et al., 2006). A part of this dust escapes the moon gravity and forms the dusty E ring of Saturn (Horányi et al., 2009;Kempf et al., 2018).
From an analysis of high phase-angle images, Porco et al. (2014) suggested a list of 100 jets of dust emission for which the coordinates and tilts were derived from images (see also Spitale et al. (2015)). To demonstrate an application of our model to the dust emission from Enceladus, we selected one single jet from this list with coordinates (−80.25 • N, 55.23 • E), which is tilted by 5 • from the surface normal in an azimuthal direction 38 • away from local north. The ejection is stationary, so that the production rate γ(t) is constant. The distributions we implemented for particle sizes, ejection speed, and direction are given by and Equation (48) was derived by Schmidt et al. (2008) to describe the acceleration of dust grains in the gas flux in the vents that supply the sources. Particles smaller than R c (measured in the same units as R) tend to accelerate up to the gas velocity (u gas ), while particles larger than R c are significantly slower. For this distribution, we have u min = 0, while u max is equal to the gas velocity u gas . Table 3 lists the parameters of the distributions and other parameters that are necessary to set up the model. Figure 12 shows the result for the number density of dust obtained from the model for the single jet, evaluated along the trajectory of Cassini during the E2 flyby. The model was evaluated for grains with a radius larger that 1.6 micron, which corresponds to the size threshold for the HRD data shown in the plot. We multiplied the model profile by a factor so that the peak matches the measured peak density. To match the HRD profile at a large distance from the plume, we added a constant background of 0.01 particles/m 3 to the model number density. The selection of the grain size in the model was realized by adjusting the function G p u appropriately (see Table 3). We obtained the position of the spacecraft from the reconstructed spice kernels of the mission (https://naif.jpl.nasa.gov/pub/naif/CASSINI/kernels/), using the NAIF Spice toolkit (https://naif.jpl.nasa.gov/naif/toolkit.html). Using these parameters, we recovered the location of the maximum number density on the Cassini trajectory (Fig. 12). Our model could now be applied to all the jets identified by (Porco et al., 2014), and the results could be fit to in situ data. Similarly, we could calculate the geometrical optical depths of the dust emitted along a given line of sight and compare this to the brightness distribution in images (see Sect. 4.3 for an example). For a quantitative comparison to images, we can apply light scattering modeling to the dust configuration that is derived from the dust distribution model. The dust ejection is stationary, the rate was chosen to fit the data Eq. (47) with µ = −1.0 and σ = 1.5 Small particles dominate the population, is used as a size distribution but the distribution is not too steep Eq. (48) with R c = 0.5 µm This expression describes and u gas = 1000 m/s is used the dust acceleration by the gas flux as the ejection speed distribution inside the channels as it is at Enceladus Eq. (49) with ω = 10 • is used The distribution describes uniform emission as the distribution into a cone of a given width of the ejection direction R min = 1.6 µm The lower sensitivity threshold of the HRD An arbitrarily chosen but reasonably high value, meaning that we can neglect particles with larger sizes. With our R max = 6 µm choice of size distribution and the size -dependent ejection speed distribution, it is highly unlikely that the HRD detects a particle of this size at the altitude of E2 flyby p = 0 We are interested in the number density

Surface deposition of material from plumes on Europa
There is evidence that cryovolcanic activity also generates plumes on the Jupiter satellite Europa Phillips et al., 2000;Roth et al., 2014;Sparks et al., 2016;Jia et al., 2018).
Owing to the higher gravity of Europa, plumes would be more confined than those on Enceladus, and they may be harder to observe or sample directly (Southworth et al., 2015;Quick and Hedman, 2020). However, surface features from plume deposits may provide evidence for past cryovolcanic activity, and if spectral features exist, allow the remote characterization of material from the interior Phillips et al., 2000;Quick and Hedman, 2020). We calculated the radial variation of the mass flux of dust grains falling back onto the surface, considering four different dust sources with identical characteristics except for the particle size distribution and the distribution of ejection directions. We considered two size distributions and two ejection modes: one describing a narrow jet, and the other describing a broader, more diffuse ejection. For the size distributions we employed a power law with two different values of the exponent q. The greater q, the more abundant the small dust particles. For the ejection direction, we used a pseudo-Gaussian distribution of the polar angle giving a nonzero probability of ejection in any direction. The normalization constant C norm was found numerically for fixed values of ψ max and ω. The dust production rate γ was identical for all the sources, meaning that they produced the same number of dust particles per unit time, but we then obtained different mass production rates for the different size distributions. For the ejection speed distribution we again used Eq. (48). Table 4 lists the distribution parameters. We can compute the rate of dust mass produced per second from the size distribution as where ρ is the density of the ice grains, and the grains were assumed to be spherical.
To compute flux instead of density, we must modify Eq. (20) as As the vector r is normal to the surface of the satellite, the factor v|cos θ| turns the density into the flux of dust falling back to the moon. Multiplying equation 53 by 4πρ/3 (ρ is the density of the particle material) and setting p = 3 and r = r M gives the mass flux onto the surface. Because r = r M on the surface, it is geometrically impossible to obtain θ * i < π/2 or any particles moving upward. The dust ejection is considered stationary Eq.( 50) with q = 3, R 1 = 0.2 µm We paired these size and ejection and R 2 = 20 µm as direction distributions the shallow size distribution to model four sources: a narrow jet and with the same boundaries, but q = 5 with a shallow size distribution, as the steep size distribution a narrow jet with a steep size distribution, a diffuse source with a shallow size Eq. (51) with ω = 5 • and ψ max = 0 • distribution, and a diffuse source for the narrow jet and with ω = 45 • with a steep size distribution and ψ max = 45 • for the diffuse source Eq. (48) with R c = 0.5 µm This expression describes and u gas = 700 m/s is used the acceleration of the dust by the gas as the ejection speed distribution flux inside the vents as it is thought be at Europa We wish to compute the total mass R min = 0.2 µm produced by each source, R max = 20 µm therefore the interval (R min , R max ) coincides with (R 1 , R 2 ) p = 3 We are interested in the mass ρ = 920 kg/m 3 Water-ice density With the parameters from Table 4, Eq (52) gives 0.6 kg of dust produced each second for the source with the shallow size distribution (q = 3) and 0.05 kg for the source with the steep size distribution (q = 5). Figure 13 shows the distribution of the dust deposition (mass flux onto the surface) with distance from the source on the moon surface. Figure 13: Dust flux on the surface of the moon at various distances from the dust source. The production rate of the sources with shallow dust size distribution is 0.6 kg/s, and the production rate of the sources with steep size distribution is 0.05 kg/s. Taking the observational evidence together, we expect the outgassing activity on Europa to be intermittent or episodic, with a potentially complex time dependence. This case could be modeled using a time-dependent function for the dust production rate. We consider a time-dependent dust ejection in our final example.

Images of a volcanic eruption on Io
The innermost satellite of Jupiter, Io, is a geologically very active body. It has multiple volcanic centers scattered over its surface (Strom et al., 1979;Keszthelyi et al., 2001;Geissler et al., 2004). The volcanic plumes vary in many features as shape, size and period of activity and often they are asymmetric. Our model allows to investigate such plumes by using inclined jets with time-dependent eruptions, or even a simultaneous dust ejection from sources with different properties.
The main flaw of our model in application to the Io volcanic plumes is that the gas emerging from the vents is generally not collisionless. In the dense cores of the plumes, condensation occurs up to a height of several hundred kilometers (Cataldo et al., 2002). Moreover, the dynamics of the dust will be influenced by the surrounding gas. Condensation (i.e., dust production) in the plume might be taken into account in principle by modeling such a plume with additional point sources placed within the condensation column, with appropriately adjusted time-dependent dust production rates. For this purpose, we implemented in the code the possibility of computing the dust density also at locations below (closer to the moon center) the location of the source. In this case, we still use the assumption that the true anomaly at the location of the source is lower than π, that is, downward ejection from a source located above the surface is not possible.
In our simple example presented here, we considered just one source located on the surface to represent a volcanic plume with a finite duration of activity. We took several images of the expanding dust plume from the same point in space with a favorable geometry when the volcano is located at the limb of the moon. We also assumed that the glow of the surface of Io was already subtracted from the images but a homogeneous background brightness was present, so that the moon disk looks dark against this background. We tilted the volcano slightly by 3 • toward local Table 5: Parameters used to construct the volcanic plume images.

Parameter
Comment Eq. (55) with t max = 500 s The volcanic plume and γ 0 = 10 14 s −1 as the time-ejected dust for 1000 s dependent dust production rate Eq. (50) with q = 3.0, R 1 = 0.2 µm, and R 2 = 20 µm is used as the size distribution The ejection speed is uniformly The narrow range of initial distributed between velocities allows us to obtain 700 m/s and 750 m/s the umbrella-shaped plume Eq. (51) with ψ max = 0 • The jets are very narrow, and ω = 5 • is used as the distribution but there is a nonzero probability of the ejection direction of ejection in any direction R min = 0.2 µm We observe the particles with the sizes R max = 0.4 µm close to the optical wavelength range, so that the radii are twice smaller p = 2 We are interested in the area covered by the particles south. The volcano was not exactly in the center of the image for the purpose of reducing the computational difficulties arising when ∆φ or ∆β ≈ 0 (see Sect. 2.4). Synthetic images of size 128x128 pixels were then constructed in the following way. Each pixel corresponds to a line of sight. We placed a grid of points on the line of sight for which we calculated the total cross section covered by the dust grains (p = 2 in Eq. (14)). Then we integrated over the lines of sight to obtain the geometrical optical depth (total particle cross section per area covered by the pixel). The color of a pixel then corresponded to the value determined for the geometric optical depth.
We employed Eq. (50) for the size distribution, Eq. (51) for the distribution of the ejection direction, and a simple size-independent expression (Eq. (54)) for a uniform ejection speed distribution, When the time interval in which the plume is observed is much shorter than the characteristic time of the plume variability, the ejection can be considered steady. However, here we modeled a nonstationary volcanic plume that was active for only 1000 seconds and constructed nine snapshots of this plume to show the evolution of the space distribution of dust with time where t max = 500 s. Table 5  value would not affect the final result. In the case considered here, the maximum ejection velocity was set to be lower than the escape velocity of Io, so that there is a maximum height that the particles can reach. Moreover, the disk of the moon covers part of the images. Taking these two facts into account, we excluded the points on the lines of sight from the calculations that crossed the moon disk and the points that lie at a greater distance from the surface than the maximum height.

Discussion
We developed a semianalytical model that allows the user to derive the average properties of dust (number or mass density, fluxes, and optical depths) ejected from an atmosphereless body. Physically, our approach is based on the two-body problem, that is, it neglects any forces on the dust particles other than the point-mass gravity of the source body. We show that the model still has a wide range of applications. These are situations where the gravity of other perturbing bodies and higher-order gravity from a nonspherical mass distribution of the central body can be neglected along the entire path that a dust particle takes from its source to its sink. The nongravitational forces must also be negligibly small, for instance, electromagnetic forces acting on charged grains, radiation-induced forces (solar radiation pressure and Poynting-Robertson drag), and drag exerted by ambient gas or plasma. For instance, the model can be applied to the Enceladus dust plume in a region that is sufficiently close to the dust sources on the south polar terrain of this satellite. It cannot be applied to estimate the dust density at higher northern latitudes of Enceladus, however, because this region is too far away from the sources near the south pole, and three-body forces due to Saturn have already affected the shape of the dust configuration. This can be seen in a peculiar pattern of color variation (Schenk et al., 2011a), that can be matched with two wedge shaped regions extending deep into the northern hemisphere (Schenk et al., 2011b;Schenk et al., 2017), for which three body models of the plume predict enhanced fall back rates of south polar dust (Kempf et al., 2010). The model can be applied to investigate the dust clouds around the Galilean moons (Krüger et al., 1999) at any longitude and latitude, however, because in this case, the dust emission occurs (nearly) uniformly over the whole surface of the satellites as long as the point of interest is located deeply enough in the Hill sphere of the satellite. Mathematically, our model relates the dust distribution at the site of ejection to the distribution at the point of interest (spacecraft position) by using the conservation laws of energy and momentum provided by the two-body problem Sremčević et al., 2003). We use the fact that the position of the spacecraft, the position of the source, and the center of the moon define the plane to which the movement is restricted. In the evaluation of the dust properties at a given point, this allows us to carry out two of the three integrations over velocity space analytically. Only one remaining integration must then be performed numerically. The distribution-based approach is very flexible, and because it allows employing asymmetric and nonstationary modes of dust emission, it allows modeling quite complex situations. A time dependence could relatively easily be introduced in the distributions of ejection speed, direction, and size as well. Furthermore, the model could also be extended to model dust emission from a comet or an active asteroid, making the Sun the central body and using a rescaled solar mass to account for the radiation pressure.
Relying only on one numerical integration, the model becomes computationally very efficient, so that even image reconstruction becomes feasible, although it involves the evaluation of the dust properties in the three-dimensional region that is in the field of view. From the examples discussed in Sect. 4, the image of the volcano (Sect. 4.3) is the most computationally demanding. Nevertheless, on a usual four-core PC, each image in Fig. 14 required 0.4 s of elapsed time to be obtained.
The model is implemented in Fortran-95, and the package called DUDI is available at https://github.com/Veyza/dudi for free usage under the terms of GNU General Public License. A user can use the probability density functions and choose from several variants that are already implemented in the described examples, or they can implement new variants.  (14) with a lognormal size distribution and a fixed value of u.

Appendix B. Physical meaning of the function G u
For each source, the first step in the numerical calculations is to obtain values of G p u (R min , R max ) (Eq. (14)) on a dense grid of u. This means that we perform integrations over the particle size R treating u as a parameter. The size distribution f R (R) can be defined on an interval of grain sizes, its lower and upper boundaries being parameters of the distribution.
These boundaries need not necessarily be equal to R min and R max . Instead, R min and R max should be understood as limits of an observable interval of particle radii. Particles smaller than R min or larger than R max may exist, and they contribute to the normalization of f R (R), but they would not contribute to the value of n (Eq. (20)). This concept is illustrated in Figs. 15 and 16. The quantity plotted in the two figures is f R (R)f u (u, R). In Fig. 15 ,f R (R) is a lognormal distribution that is formally defined in the interval (0, +∞), or in other words, +∞ 0 f R (R)dR = 1, and the interval (R min , R max ) does not cover the whole domain of f R (R). In Fig. 16, the particle sizes are distributed between certain R 1 and R 2 as Eq. (50). In this example, R min = R 1 , but R max < R 2 . If R min were lower than R 1 , it would not have changed the result because for R < R 1 , we have f R (R) = 0.
Although we defined the size distribution outside the interval (R min , R max ), its shape there does not affect the final result as long as it remains the same inside (R min , R max ). However, we can investigate the defined size distribution by applying different values of R min and R max . This approach implies that the size distribution and the initial speed distribution are physical characteristics of the dust source, while R min and R max represent the sensitivity range of the instrument with which we perform observations. In the model implementation DUDI, the interval (R min , R max ) can coincide with the f R (R) domain or can be even larger. The units of the particle radius R matter only in the definition of f R and f u , so that we suggest measuring R in microns to have simpler numbers in the expressions. With different formulae for f R and f u , shorter expressions for G p u can be obtained that require less cumbersome computations (see, e.g., Postberg et al. (2011)). However, we purposefully consider a set of f R and f u in our model that cannot be simplified to show the general form of the solution and to allow flexibility in applying the model.