Peripheral Collisions of Ice-covered Silica Dust Grains

Collisions with ice-covered silica grains are studied using molecular-dynamics simulation, with a focus on the influence of the impact parameter on the collision dynamics. The ice mantle induces an attractive interaction between the colliding grains, which is caused by the melting of the mantles in the collision zone and their fusion. For noncentral collisions, this attractive interaction leads to a deflection of the grain trajectories and, at smaller velocities, to the agglomeration (“sticking”) of the colliding grains. The bouncing velocity, which is defined as the smallest velocity at which grains bounce off each other rather than stick, shows only a negligible dependence on the impact parameter. Close to the bouncing velocity, a temporary bridge builds up between the colliding grains, which, however, ruptures when the collided grains separate and relaxes to the grains. At higher velocities, the ice in the collision zone is squeezed out from between the silica cores, forming an expanding disk, which ultimately tears and dissolves into a multitude of small droplets. An essential fraction of the ice cover in the collision zone is then set free to space. Astrophysical implications include the possibility that organic species that might be present in small concentrations on the ice surface or at the ice–silica interface are liberated to space in such noncentral collisions.


Introduction
In interstellar and interplanetary environments, nanoparticles ("dust grains") are abundant. Protoplanetary dust disks (Armitage 2010(Armitage , 2011Blum 2010;Birnstiel et al. 2016), debris disks of developed planetary systems (Nakamura et al. 2012;Gáspár et al. 2013), planetary rings (Esposito 2010), and the dust tails of comets (Bentley et al. 2016;Langevin et al. 2016) provide well-known examples. Beyond the ice line, not only pure ice particles but also inhomogeneous grains are found, consisting of a silica core covered by an ice mantle. The population of dust grains changes due to collisions; if colliding grains stick together, they form aggregates, which form the building material for the further evolution of the dust cloud. At larger velocities, grains bounce off each other or may even fragment. The threshold velocity-called "bouncing velocity"-is thus an important characteristic of dust grains (Bridges et al. 1996;Dominik & Tielens 1997).
The bouncing velocity strongly depends on the material from which grains are built. Available macroscopic formulae for homogeneous grains (Johnson et al. 1971; Thornton & Ning 1998;Brilliantov et al. 2007;Krijt et al. 2013) predict it to increase with surface energy and to decrease with the Young modulus, the mass density, and the radius; these findings are in agreement with experiment (Bridges et al. 1996;Higa et al. 1998;Schäfer et al. 2007;Krijt et al. 2013;Heim et al. 1999;Poppe et al. 2000;Güttler et al. 2013). Atomistic simulations of nanoparticle collisions also show these dependencies for, among others, silica grains (Nietiadi et al. 2017a(Nietiadi et al. , 2019; in addition, they show that phase transformations in the collision zone may alter the behavior (Nietiadi et al. 2017a(Nietiadi et al. , 2020b. Thus, water ice may melt and sinter the colliding grains together (Nietiadi et al. 2017a).
The effects of the ice mantle surrounding dust grains on dust aggregation may be pronounced. Kimura et al. (2020) argue that beyond the snow line, the ice mantle helps aggregation, and sublimated water ice inside the snow line may serve a similar purpose. Wang et al. (2005) reason that an ice coating may hold porous grains together in collisions. Korablev et al. (2021) recently proposed that salt stored in dust grains may interact with water and radiation to generate HCl, thus pointing at the relevance of water and ice mantles in planetary environments.
The collision behavior of ice-mantled silica grains is particularly intriguing as no macroscopic description exists for such inhomogeneous collisions. Previous simulations studied only central collisions and revealed the influence of mantle melting on grain sticking (Nietiadi et al. 2020a). However, central collisions are the exception, as usually the impact parameter will not vanish; peripheral collisions are the rule. In the present study, we explore to what extent an ice mantle influences off-center collisions of silica-core ice-mantle grains.

Method
The core-mantle grains studied here consist of a silica core of radius R = 20 nm covered by a water-ice mantle of thickness d = 5 nm. Both the core and the mantle are amorphous. The silica interaction is described by the Munetoh et al. (2007) potential, while the water ice is simulated using the Molinero & Moore (2009) potential. The interaction between silica and water is modeled by Lennard-Jones potentials; details on the interaction potentials and the construction of the core-mantle grains are provided in Nietiadi et al. (2020a). Initially, the grains are at a temperature of 150 K. The molecular-dynamics simulations are performed with the LAMMPS code (Plimpton 1995). Atomistic snapshots are generated with OVITO (Stukowski 2010).

Collision Dynamics
In the following, we will analyze the collision processes for two impact parameters, b = 40 nm and b = 20 nm, which we will denote as glancing and oblique collisions, respectively. We consider impact speeds of 100-1000 m s −1 . These speeds are much larger than those thought to take place in Saturn's rings (Hatzes et al. 1988) and protoplanetary disks, especially for 25 nm grains (see Ormel & Cuzzi 2007;Windmark et al. 2012). However, such high speeds might occur in debris disks of developed planetary systems (Krivov et al. 2006;Gáspár et al. 2012). Figure 1 gives an example of the collision dynamics for a glancing collision at an impact parameter of b = 40 nm; at this impact parameter, the naked silica cores would just miss each other, such that the interaction between the ice mantles governs the collision dynamics. The velocity is just above the bouncing threshold. At the moment of closest approach, Figure 1(b), the ice mantles overlap, while the silica cores are clearly separated; this means that the two grains repelled each other; their distance is 43.95 nm at the closest approach, 76 ps. After passing each other, a water bridge joins the two grains, Figures 1(d) and (e), but breaks eventually, Figure 1(f). At velocities of 165 m s −1 and below, the two grains remain bound together and do not separate, but circle each other. We therefore determine the bouncing velocity as 170 m s −1 for this case; it was determined with a 5 m s −1 accuracy.
The case of an impact parameter of b = 20 nm is illustrated in Figure 2. At this impact parameter, the two silica cores would touch directly, if they were not shielded by the ice mantle. During maximum compression, Figure 2(b), corresponding to the closest distance (38.21 nm) of the two grains, the silica cores are still covered by a water mantle even in the contact zone; but due to the high pressures occurring here, the ice is squeezed out laterally. This phenomenon is hardly observable for the glancing collision, b = 40 nm. The two grains then repel each other; upon separation, a wide water bridge is established between them, which has a hollow interior; see Figures 2(c)-(e). The bridge finally breaks, Figure 2(e). Note that the final velocity vectors show that in the late part of the interaction, the water bridge induced an attractive interaction between the grains and reduced the deflection angle.
Quantitative information about the collision dynamics is obtained by measuring the coefficient of restitution (COR) and the deflection angle. The COR is determined from the relative velocity of the grains after the collision, v¢, and the relative velocity before the collision, v, as Sticking collisions are therefore characterized by e = 0. Analogously, the deflection angle, θ, is the angle between v¢ and v; we give it a positive sign for repulsive interaction, i.e., if in snapshots like Figure 2, the grain approaching from the right-hand side is deflected upwards, and the grain approaching from the left-hand side is deflected downwards. For sticking collisions, θ is undefined. For the glancing collision, the deflection angle slowly converges to 0 at high velocities, Figure 3(a). For smaller velocities, the deflection angle changes sign; at v = 220 m s −1 , it is zero. This change means that apart from the repulsive interaction dominating at large velocities, at small velocities, an attractive interaction also sets in, which leads to negative deflection angles. This attractive interaction originates from the fusion of the ice mantles of the two grains. At the smallest bouncing velocity monitored by us, 170 m s −1 , it is e = 0.60 and θ = −37.2°.
The deflection angle for the oblique collision, Figure 3(b), has a more pronounced velocity dependence than for glancing angles in that the maximum repulsive deflection angle of 88.8°a ttained at a velocity of v = 400 m s −1 is considerably larger. At higher velocities, grain repulsion becomes less effective due to the high collision velocity. At smaller velocities, the attractive action of the ice mantle and the formation of the connecting water bridge lower the deflection angle again. The collision shown in Figure 2 has a deflection angle of 26.4°. We plot in Figure 4 the velocity dependence of the COR for the cases studied here and compare it to that for central impact (Nietiadi et al. 2020a). The results can be rationalized by comparing them to the law ), which has been derived for collisions of adhesive spheres within the framework of the Johnson-Kendall-Roberts theory of adhesive contacts (Johnson et al. 1971;Thornton & Ning 1998;Brilliantov et al. 2007). v b denotes the bouncing velocity, and the prefactor α takes a reduction of the COR by dissipative processes into account. Such dissipative processes include collision-induced heating and excitation of grain vibrations for homogeneous grains (Chokshi et al. 1993;Dominik & Tielens 1997;Krijt et al. 2013); in our case, plastic deformation and phase transformation of the ice mantle also contribute (Nietiadi et al. 2017a(Nietiadi et al. , 2020a. Figure 4 shows the evolution of the COR with impact velocity v for glancing collisions (b = 40 nm). For high velocities, the COR approaches 1, as energy losses in this glancing collision become minor. In this case, the fit, Equation (2), works fine with α = 1; only at small velocities around the bouncing velocity do the simulation data show a more abrupt decline of the COR than the theoretical law, Equation (2). For smaller impact parameters, α has to be reduced to 0.75 and 0.65 for b = 20 nm and 0,   respectively, signaling increased energy losses during the more central collisions. While the simulation data for the oblique collision with b = 20 nm follow the law, Equation (2), quite satisfactorily, the data for central impact show some deviations.
Here the COR exhibits a broad maximum of e = 0.63 at around 600 m s −1 ; toward higher velocities, the COR slightly decreases again. Note that at this velocity, for b = 20 nm, the attractive interaction of the ice mantles is about to cancel the repulsion of the silica cores. We indeed find that for velocities above 400 m s −1 , the silica cores are increasingly compressed and thus contribute to energy dissipation.
The fit of the COR to Equation (2) allows us to obtain the bouncing velocity with higher precision than the 5 m s −1 interval of our simulation results. It gives us v b = 195.0 (203.8, 164.1) m s −1 for b = 0 (20, 40) nm. The variation of the bouncing velocity with impact parameter is thus astonishingly small, only 16% even for glancing collisions.
This may be compared with collisions of naked silica spheres. Here the bouncing collision for central impacts amounts to 469 m s −1 (Nietiadi et al. 2017b). Interestingly, this does not change much for oblique collisions. For b = 20 nm, we find that the grains stick at 470 m s −1 and bounce at 475 m s −1 . We conclude that the high adhesion of naked silica leads to high bouncing velocities with little dependence on the impact parameter.
This situation is in contrast to Lennard-Jones clusters, where a strong dependence of v b on b has been found (Kalweit & Drikakis 2006). Here, v b decreases by more than a factor of 2 if b increases from 0 to R/2 (corresponding to b = 20 nm in our case); for more glancing collisions, v b decreases further to zero.

Deformation of the Ice Mantle and Water Ejection
The silica cores deform only negligibly during the collision at small velocities-see Figures 1 and 2-and even at higher velocities the deformations are elastic and relax after the collision. In contrast, the ice mantles are strongly affected by the collision, and the ice exhibits high plasticity. Already in the vicinity of the bouncing velocity, the high pressure acting during the phase of the closest approach squeezes out the ice laterally; the induced bulge in the collision plane is particularly well visible for b = 20 nm in Figure 2(b) but relaxes later on.
This high plasticity of the ice mantle is caused by the collision-induced heating of the ice and by the high shear stresses acting. We determined the temperature distribution in a 10 Å wide slab around the collision plane for the bouncing collisions depicted in Figures 1 and 2 and found that at the moment of closest approach, 28% of the water molecules are above the triple-point temperature for b = 20 nm, and 7% for the less violent collision for b = 40 nm. Thus, even though the collision zone is not entirely molten, the material is considerably easier to deform than at the initial temperature of 150 K.
At the closest approach, strong shear builds up in the grains. For the oblique collision with b = 20 nm, the shear strain is displayed in Figure 5(a). Shear is largest at the interface, as expected. Only a thin region of the silica core is strongly affected, while a large volume of water is displaced and displays high shear. The hemispherical shear pattern in the silica core is similar to the one expected from contact in isotropic media, with a maximum below the contact area (Johnson 1985). The volumetric strain in the contact area is shown in Figure 5(b). Large values-reaching strains of 1 and larger-are dominant in the central region of the ice mantle where the ice has been squeezed out laterally. Values in silica, however, are moderate; in particular, they are below the threshold of 20%, above which an amorphous-amorphous phase transition may occur (Dávila et al. 2003). Note that in neither Figures 5(a) nor (b) are shear bands observed in silica; this signals that the deformation in the silica core is elastic and the silica will return to its initial structure after the collision.
The glancing collision at an impact parameter of b = 40 nm produces almost no appreciable shear in the silica cores, Figure 5(c), except an extremely thin region next to the water. The ice mantle undergoes considerable shear strain in the channel between the colliding grains, as it is squeezed out. The shear in the ice is diminished in the vicinity of the ice-core interface.
At small velocities above the bouncing velocity, the ice forms an adhesive neck between the colliding grains. During the outgoing trajectories of the collided grains, this neck becomes elongated and hence decreases its thickness. This process is particularly well seen in snapshots (c)-(e) of Figure 2, which show that the initially hollow neck collapses laterally. Finally, the bridge formed between the two departing grains breaks, Figures 1(f) and 2(f), and its remainders of the bridge are absorbed again in the ice mantles of the nanoparticles. Note that in these two events, no water molecule was ejected; i.e., the rupture of the bridge occurred sufficiently softly so that all water molecules eventually were reintroduced into the ice mantle. The ice mantles have restored their original shape after the collision. At higher collision velocities, though, another process takes over. The temporary bump seen in the collision zone becomes more energetic and evolves into a disk-like structure that expands from the collided grains. This process is displayed in Figure 6 for the oblique collision, b = 20 nm. Figure 6(a) shows the nearly symmetric expansion of the disk. Later, it shatters into a ring that is still connected by several spokes to the two grains, Figure 6(b). These spokes fragment earliest, while the ring survives in the final snapshot that we show, Figure 6(c). It already starts decaying to a series of droplets by Rayleigh instability; this process, however, will take more time. Figure 6(b) in addition shows a bridge directly connecting the two grains, which is ruptured in Figure 6(c); however, in contrast to lower velocities, this bridge now contains only a small fraction of the ejecta.
For b = 20 nm, the formation of a collision disk starts at 500 m s −1 and governs the evolution for all larger velocities. At 500 m s −1 , Figure 6, around 210,000 water molecules are ejected from the grains-most of them are contained in the disk structure developing into ejected droplets. The number of ejecta continuously increases with collision velocity until they reach ∼400,000 water molecules at v = 1000 m s −1 . This number constitutes the majority of water contained in the collision zone, as can be seen as follows. The collision zone can be approximated by a cylinder of height 2d and radius R; see Figure 2(b). With the molecular density of water of n = 33.43 nm −3 , the number of water molecules contained in the collision zone is n × 2πdR 2 , amounting to 420,000 molecules. This is a small fraction of the total number of water molecules in the two grains, N tot = 2.25 × 10 6 .
For glancing collisions, such a disk structure also develops; however, due to the small deflection angles, Figure 3, the two collided grains remain in the plane in which the disk expands. Figure 7(a) shows the asymmetric expansion of the disk caused by the two grains dragging the disk in opposite directions within its plane. While the disk is not yet ruptured in Figure 7(b), it has torn into a series of bridges connecting the two departing grains in Figure 7(c). In the further evolution of this structure, the bridges will rupture and develop into a series of ice droplets. Note that the collided surfaces look strongly altered by the collision, in contrast to the case of oblique collision at b = 20 nm, Figure 6.
For the glancing impact parameter, we observe at v = 500 m s −1 still the development of a single bridge between the departing grains. Disk-like structures start to develop around 700 m s −1 and become dominant at the case shown in Figure 7, 900 m s −1 . The number of ejected molecules at the end of the simulation is still rather modest, amounting to only 7000 in the case of Figure 7. This will, however, change when continuing the simulation until the bridges tear, and numbers similar to those in the case of the oblique collision, b = 20 nm, are to be expected. Figure 8 summarizes the evolution of the number of ejected water molecules with collision velocity. As discussed above, oblique collisions lead to considerably more intense ejection activity than glancing collisions. The onset of ejecta formation is at around 500 m s −1 and coincides with the formation of a collision disk for b = 20 nm, while for b = 40 nm, until the end of the simulation, ejection activity does not feature a pronounced velocity dependence. The silica cores remain unaffected by the collision at the velocities studied here, v 1 km s −1 because the cohesive energy of silica is considerably larger than that of water ice. This result applies to the destruction of single grains by high-speed collisions; it differs from the assumptions of much lower grain fragmentation speed in the modeling of dust-aggregate evolution in protoplanetary disks (Brauer et al. 2008;Birnstiel et al. 2010a); see our discussion in Section 3.3.
The angular momentum carried away by the ejecta is modest, in agreement with the small fraction of mass they  carry. In the case of maximum ejection observed in the present study-oblique collision at v = 1 km s −1 ; see Figure 8, where 6.5% of the water molecules were ejected-the ejecta carry 3.5% of the original angular momentum. In all other cases, the angular momentum loss is smaller.

Discussion
It may appear astonishing that the ice mantle does not fracture during this violent collision. Besides the collision-induced heating effect mentioned above, a size effect may also be relevant. Many brittle materials show ductile behavior at the nanoscale. For silica, this effect is well known and has been related to an increase in the rate of bond-switching events near surfaces, leading to a high fracture strength (Luo et al. 2016); this effect is emphasized at high strain rates (Ramachandramoorthy et al. 2019). We assume that a similar mechanism may be operative in amorphous water ice to render it more ductile at the nanoscale. Yasui et al. (2017) and Wong & Baud (2012) discuss experimental work on a related effect-the ductile-to-brittle transition in ice-silica mixtures as a function of composition and in porous rocks-but on a macroscopic scale.
The immense water ejection yields in noncentral collisions may be important for astrophysical implications. It has been assumed that the surface of water ice in a space environment may host organic molecules, which may be altered by UV irradiation from nearby stars and cosmic rays (Ehrenfreund & Charnley 2000;Tielens 2005;Bennett et al. 2013;Gudipati & Castillo-Rogez 2013;Cottin et al. 2017). Also, the silica-ice interface could be a host of organics; this interface is special because energy deposited in the grain by cosmic-ray irradiation will primarily heat the silica-due to its larger stopping power -such that the liberated heat may induce reactions of organic molecules or radicals in the vicinity of the interface (Bringa & Johnson 2004;Arumainayagam et al. 2019). These organic molecules-both at the surface and near the interface-may be freed from the grains by the oblique collisions described here to be ejected into space and be detected by IR astronomy (Johnson 1990;Li & Draine 2001;Bringa & Johnson 2002). As long as such organic molecules are present only in small concentrations, they will not change the collision dynamics, which is governed by the silica core and the ice mantle; if they present a majority species, their influence on the collision dynamics may be more profound.
This process may be considerably more efficient than sputter emission by ion irradiation from cosmic rays or other incident particle radiation (Barlow 1978a(Barlow , 1978bVasylyev et al. 2004).
Here, it has been shown that such sputter events preferably lead to the fragmentation of the ejecta because the imparted energies lead to bond breaking in the organics (Kitazoe & Yamamura 1980;Bitensky & Parilis 1987;Urbassek 1987;Anders et al. 2020). In contrast, the emission by glancing collisions is sufficiently soft to guarantee the ejection of the intact molecules.
The dynamics of colliding ice-covered silica grains is fundamentally different from that of pure ice nanograins, as these latter never show bouncing (Nietiadi et al. 2017a) because the two nanoparticles (NPs) melt and sinter together. Pure ice grains may be considered as the limiting case of large mantle thickness, d ? R. In this sense, the present results on ice-mantled grains interpolate between collisions of pure silica (d = 0 nm) and pure ice NPs (R = 0 nm).
Our grains were initially at a temperature of 150 K, close to the temperature of the snow line of water ice (Sasselov & Lecar 2000;Podolak & Zucker 2004). As temperature decreases with an inverse power law with distance from the star in protoplanetary disks (Andrews & Williams 2007;Öberg et al. 2011)-and also in our solar system (Lewis 1974)-this condition is only valid at a certain distance from the star. Water ice, if condensed or adsorbed in astrophysical environments at such low temperatures, will be amorphous (Baragiola 2003;Baragiola et al. 2013;Gärtner et al. 2017).
The initial grain temperature will influence the collision outcome in several ways. The most immediate influence might be that a higher initial temperature requires less collisional heating-and thus smaller collision velocities-to bring the ice mantles above the melting temperature, resulting in a softer and more sticking contact. In addition, the elastic modulus and the specific surface energy of ice will decrease for higher initial temperature, again slightly influencing the collision dynamics.
The influence of the mantle thickness d on the collision dynamics of core-mantle grains was studied previously for central collisions (Nietiadi et al. 2020a). There, d was systematically varied, and it was shown that a larger mantle thickness increases the bouncing velocity and decreases the COR for bouncing collisions. Both effects originate from the "cushion"like damping of the ice mantle that dissipates energy and increases adhesion (sticking). It may be expected that these effects will also carry over to the case of noncentral collisions.
For central collisions, the effect of the core size, R, on the collision dynamics was also investigated (Nietiadi et al. 2020a). In agreement with available macroscopic collision theories (Johnson et al. 1971), increasing R leads to a decrease in the bouncing velocity. This feature holds true even though collision theories are designed for homogeneous grains, and the simulations are for inhomogeneous core-mantle grains. The effect of the adhesive mantle can be thought of as increased effective surface energy, which increases with mantle thickness d.
In the present study, we chose an idealized model scenario for our simulations: a spherical silica core is surrounded by an ice mantle of constant thickness. Reality may deviate in several respects from this model scenario: the core and/or mantle may not be spherical, the mantle thickness may vary with position, core and mantle material may be mixed together, etc. Such deviations may occur during the formation process of the grains, but may also be altered by previous collisions that the grains experienced in their history. In this sense, the simulations reported here are meant to illustrate the principal processes under which an ice mantle modifies the collision of bare silica cores. At velocities higher than those considered in the present study, i.e., above 1 km s −1 , grain fragmentation processes will become dominant. These already showed up in the disk-like structures originating from the ice mantles that were studied in Section 3.2, but will also involve the silica cores at higher velocities and ultimately destroy the grains. In the extreme limit when the collision velocity exceeds 1 km s −1 , shock-wave grain densification will be important (Kubota et al. 2002;Barmes et al. 2006), and for higher velocities, leading to pressures above ∼35 GPa, phase changes might occur within the silica grains (Tracy et al. 2018).
Finally, discussions of dust collisions in an astrophysical context are often based on collisions of grain agglomerates (Dominik & Tielens 1997;Krause & Blum 2004;Paszun & Dominik 2006;Blum & Wurm 2008); i.e., each collision partner is itself not a simple grain but has a complex structure composed of an assembly of individual grains held together by surface forces (Blum 2010). Such aggregates are thought to possess a high porosity because each individual grain is in contact with only a few neighboring grains. The collision dynamics of such aggregates deviates strongly from the collision dynamics of individual grains because energy dissipation may not only occur within the grains but in the contacts between grains. Thus, at low velocities, three orders of magnitude lower than those simulated here (Blum & Wurm 2008), aggregate compression and restructuring are typical collision results; also, aggregate fragmentation is considerably easier to achieve than in grain collisions because the surface forces between grains are considerably smaller than the surface forces between atoms (Dominik & Tielens 1997;Ringl et al. 2012;Planes et al. 2021).

Conclusions
The collision dynamics of ice-covered silica grains deviates strongly from that of bare silica grains. The ice mantle induces an attractive interaction between the colliding grains, which influences the bouncing velocity of the grains. It also alters the collision dynamics by causing an attractive deflection of the grain trajectories, which leads-for small collision velocitiesto grain sticking. The bouncing velocity is found to be rather independent of the impact parameter.
The ice mantle exhibits enhanced plasticity caused by collision-induced heating and high strains acting. Close to the bouncing velocity, a temporary bridge builds up between the colliding grains, which, however, ruptures when the collided grains separate and relaxes to the grains. At higher velocities, the ice in the collision zone is squeezed out from between the silica cores, forming an expanding disk, which ultimately tears and dissolves into a multitude of small droplets. A large fraction of the ice cover in the collision zone is then set free to space. Astrophysical implications include the possibility that organic species that have formed in small concentrations on the ice surface or at the ice-silica interface are liberated to space in these noncentral collisions.
In addition, our results demonstrate that core-mantle grains and in particular noncentral collisions open up new scenarios that may differ in several respects from the central collisions of homogeneous grains. While the sticking behavior may be similar, bouncing collisions follow another (attractive) dynamics and may eject a considerable amount of matter to space. These features may be relevant for theoretical approaches that model the evolution of dust clouds through particle collisions (Birnstiel et al. 2010a(Birnstiel et al. , 2010b), e.g., by granular mechanics codes (Dominik & Tielens 1997;Suyama et al. 2008;Wada et al. 2008;Ormel et al. 2009;Wada et al. 2011) or by Monte Carlo algorithms (Zsom et al. 2011;Birnstiel et al. 2016).