Dynamics of irregularly shaped cometary particles subjected to outflowing gas and solar radiative forces and torques

MNRAS 510, 5142–5153 (2022) https://doi.org/10.1093/mnras/stab3769 Advance Access publication 2021 December 27 Dynamics of irregularly shaped cometary particles subjected to outflowing gas and solar radiati v e forces and torques Fernando Moreno , 1 ‹Daniel Guirado , 1 ‹Olga Mu ˜ noz, 1 ‹Vladimir Zakharov, 2 Stavro Ivanovski, 3 Marco Fulle , 3 Alessandra Rotundi, 2 Elisa Frattin 4 and Ivano Bertini 4 1 Instituto de Astrofisica de Andalucia, CSIC, Glorieta de la Astronomia s/n, E-18008 Granada, Spain 2 INAF - Istituto di Astrofisica e Planetologia Spaziali, Ricerca Tor Vergata, Via Fosso del Cavaliere 100, I-00133 Rome, Italy 3 INAF - Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34143 Trieste, Italy 4 Dipartimento di Astronomia, Universita’ di Padova, vicolo Osservatorio 2, I-35122 Padova, Italy Accepted 2021 December 22. Received 2021 December 22; in original form 2021 October 1 A B S T R A C T The dynamics of irregularly shaped particles subjected to the combined effect of gas drag and radiative forces and torques in a cometary environment is investigated. The equations of motion are integrated over distances from the nucleus surface up to distances where the gas drag is negligible. The aerodynamic forces and torques are computed assuming a spherically symmetric expanding gas. The calculations are limited to particle sizes in the geometric optics limit, which is the range of validity of our radiative torque calculations. The dynamical behaviour of irregular particles is quite different to those exhibited by non-spherical but symmetric particles such as spheroids. An application of the dynamical model to comet 67P/Churyumov–Gerasimenko, the target of the Rosetta mission, is made. We found that, for particle sizes larger than ∼10 μm, the radiative torques are negligible in comparison with the gas-driven torques up to a distance of ∼100 km from the nucleus surface. The rotation frequencies of the particles depend on their size, shape, and the heliocentric distance, while the terminal velocities, being also dependent on size and heliocentric distance, show only a very weak dependence on particle shape. The ratio of the sum of the particles projected areas in the sun-to-comet direction to that of the sum of the particles projected areas in any direction perpendicular to it is nearly unity, indicating that the interpretation of the observed u-shaped scattering phase function by Rosetta /OSIRIS on comet 67P coma cannot be linked to mechanical alignment of the particles. Key words: methods: numerical – comets: individual: 67P. 1 INTRODUCTION In previous papers, (Ivanovski et al. 2017a , b ) have provided a detailed analysis of the dynamics of spheroidal dust particles in the vicinity of a cometary nucleus by assuming a gas model characterised by a spherically symmetric expanding flo w. Iv anovski et al. ( 2017a ) pro v ed that the dynamics of such aspherical particles is markedly different to spherical particles of the same volume equivalent radius. The difference in particle shape and initial orientation on the nucleus surface leads to velocity dispersion, and the maximum liftable size is minimum for spherical particles with respect to spheroidal ones. The spheroidal particles, after some time since ejection from the nucleus, called t rot , start to experience full rotation, which depends on the particle physical parameters and nucleus outgassing properties. The model has been applied to the analysis of Rosetta /GIADA data on comet 67P/Churyumov–Gerasimenko (Ivanovski et al. 2017b ). They found that the GIADA data are best reproduced with oblate particles rather than prolate. On the other hand, Fulle et al. ( 2015 ), using the particle tracks as imaged by Rosetta /OSIRIS and a non-spherical  E-mail: [email protected] (FM); [email protected] (DG); [email protected] (OM) particle model, found that the best agreement between measured and computed frequencies was reached for oblate spheroids as well. Based on the argument that the gas drag and the gravitational attraction of the nucleus constitute the dominating forces at distances close to the comet nucleus, Ivanovski et al. ( 2017a ) neglected the solar radiation pressure force and torque in their model. Ho we ver, the effect of the solar radiative force and torque need to be estimated as they might become important at larger ( R  5 R N , where R N is the nuclear radius) distances. We hav e dev eloped a model including the gas drag, the gravitational attraction of the nucleus, and the solar radiation pressure to asses the long-term evolution of the particles under the combined effect of such forces. Instead of spheroidal particles, we have considered more realistic irregularly shaped particles having a variety of axes ratios. Irregularly shaped particles have also been used by ˇ Capek ( 2014 ) to study meteoroid dynamics. Our model is based on the rigorous integration of the equation of motion simultaneously with the Euler dynamical equations, including the effects of both gas and radiation forces and torques. The gravitational torques are not included, as they have been shown to be negligible by Ivanovski et al. ( 2017a ). In Section 2, we provide a detailed description of the equations involved in the aerodynamic and radiative force and torque models, and ©The Author(s) 2021. Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://cr eativecommons.or g/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 Dynamics of cometary irregular particles 5143 the dynamical equations, where numerical integration is performed using a quaternion-based scheme. A validation of the computer code is made through comparison with the results obtained by Ivanovski et al. ( 2017a ), by switching-off the solar radiative force and torque. In Section 3, results of the full model including aerodynamic and radiative forces and torques are provided, as well as a discussion on the combined effect of such forces. We focus on the evolution of frequencies, degree of tumbling, direction of the angular momentum v ectors, and v elocities. We also inv estigate the feasibility of the hypothesis of alignment of the largest particle surface areas respect to the solar radiation, which would lead to backscattering enhancement. Finally, Section 4 lists the conclusions that can be drawn from the present study. 2 THE DYNAMICAL MODEL This section is divided into three parts. The first part describes the dust particle models, the second one the outflowing gas model, and the third part, the radiative force and torque model. As indicated abo v e, no additional effects, such as the gravitational torque, are taking into account. 2.1 Dust particle model build-up Irregular shape model particles were adopted from the 3D Asteroid Catalogue ( ht tps://3d-ast eroids.space/), which contains 3D models of known minor bodies derived from light-curve inversion. We searched the data base to find three distinct kind of shapes, namely flattened, elongated, and rounded, to explore the differences in their dynamical behaviour. We selected the shape models of asteroids (943) Begonia , (857) Glasenappia , and (94) Aurora as representative of those dif ferent shapes, respecti vely. Vie ws of those shape models displayed with MeshLab (Cignoni et al. 2008 ) are given in Fig. 1 . In addition, spheroidal model particles have also been generated in order to perform the code validation through comparison with the results by Ivanovski et al. ( 2017a ) (see Subsection 2.4.3). In all cases, the surfaces are build-up by triangular meshes. We assume that the particles behave as solid rigid bodies, they are homogeneous, and do not contain any volatiles, so that no rocket forces are applied, and are isothermal, being characterised by a constant temperature T d . The particle density is set to ρd = 800 kg m −3 , which is consistent with OSIRIS and GIADA estimates for comet 67P (Fulle et al. 2016 ). The particle mass is denoted by m d , related to its ef fecti ve radius, r eff by m d = 4 3 πρd r 3 eff . 2.2 Aer odynamic for ce and tor que For the gas model, we followed the model by Ivanovski et al. ( 2017a ). Briefly, a non-rotating, spherical nucleus, of radius R n and mass M n is located at a certain heliocentric distance r h . The nucleus surface temperature is denoted by T s , and is emitting water molecules at a production rate of Q g . The gas is assumed to behave as an ideal gas expanding into vacuum, so that it has initial sonic v elocity giv en by V s =  γT s k B /m g where γis the ratio of specific heats, m g is the mass of a water molecule, and k B is the Boltzmann constant. The gas density ( ρg ), velocity ( V g ), and temperature ( T g ) are computed as a function of the distance to the nucleus, r , by using the analytical e xpressions giv en by Zakharo v et al. ( 2018 ), appropriate for an adiabatic spherical expansion. The aerodynamic force on the particle is given by (see Ivanovski et al. 2017a ): F a = − S ( p ˆ n + τ( V  r ׈ n ) ׈ n / | ( V r ׈ n ) ׈ n | ) d S, (1) where p and τare the pressure and the shear stress of the surface element d S , ˆ n is a unit vector along the outer normal to the surface element, and V r is the gas to dust relative velocity. As already shown by Ivanovski et al. ( 2017a ), under the typical physical parameters assumed for the nucleus and the particles, the flow o v er the particles can be considered as free molecular, and the mean collision rate of gas molecules with a dust particle is al w ays much higher than the rotation frequency of the particles. The free molecular expressions for p and τcan be found in e.g. Bird ( 1994 ), and are given here for completeness as: p/p g =  s  cos β/ √ π+ 1 2  T d T g  exp ( −s  2 cos 2 β) +  1 2 + s  2 cos 2 β+ 1 2 √ πs  cos β T d T g  ×1 + erf ( s  cos β) (2) τ/p g = s  sin β/ √ π × exp ( −s  2 cos 2 β) + √ πs  cos β1 + erf ( s  cos β)  . (3) In these equations, p g = ρV r  2 /(2 s  2 ), s  = V r   m/ (2 k B T g ) , and V r  = V g −V d  , where V d  is the velocity of the surface element taking into account the rotation of the particle, and βis the angle between V r  and −ˆ n . Finally, the aerodynamic torque is given by: M a = − S l ×( p ˆ n + τ( V  r ׈ n ) ׈ n / | ( V r ׈ n ) ׈ n | ) d S, (4) where l is the vector from the particle centre of mass to the surface element d S . 2.3 Radiati v e force and torque The calculation of radiative forces and torques on non-spherical particles is a very CPU and memory consuming task. These quantities can be computed using light-scattering codes such as the Discrete Dipole Approximation (the DDA code, see Draine & Flatau 1994 ). The problem arises when particles of large size parameter are considered, as the DDA method requires a very large number of dipoles to build up the scatterer (consequently a large memory), and needs a very large CPU time even for a single orientation of the particle, so that in practice it becomes inefficient to describe the dynamics of a given particle for large integration times. To compute the radiative force and torque on the particles, knowing that their dimensions are al w ays much larger than the wavelength of the incident light (assumed at λ= 0.6 μm), we applied a geometric optics approximation. We adopted the expressions used by Beletskii ( 1966 ) to compute the force and torque on artificial satellites by the solar radiation, as follows. The radiation pressure at a distance r h from the Sun is given by p r = ¯ k r 2 h , where ¯ k = E s 4 πc = 1 . 01 ×10 17 kg m s −2 , being E s = 3.8 ×10 26 W, the total power radiated by the Sun, and c the speed of light. The solar radiation force on the particle is MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 5144 F. Moreno et al. Figure 1. Model particles used in the dynamical calculations taken from the 3D Asteroid Catalogue ( ht tps://3d-ast eroids.space/). On the top row of panels, (943) Begonia shape model particle, having a flattened shape, on the middle row (857) Glasenappia shape model particle (elongated), and on the bottom row, (94) Aurora shape model particle displaying a rounded shape. Approximate outer dimensions of the particles relative to x -axis are (1.00 ×0.72 ×0.36), (1.00 ×0.43 ×0.36), and (1.00 ×1.09 ×1.14), respectively. given by: F r = −p r (1 −0 ) ˆ t  S ( ˆ t ·ˆ n )d S −2 0 p r  S ˆ n ( ˆ t ·ˆ n ) 2 d S, (5) where ˆ n is a unit vector along the outer normal to the surface element, ˆ t is a unit vector in the opposite direction to the solar flux, d S is an elementary surface area on the illuminated fraction of the particle, and 0 is the reflection coefficient. The first term in equation (5) corresponds to the force e x erted by the incident flux and the second term is the force produced by the reflected flux. The reflection coefficient is computed through the Fresnel equations for a given comple x refractiv e inde x m = n r + in i . Since the incident solar light is unpolarized, the reflection coefficient is calculated as the average of the coefficients corresponding to the s and p polarizations (see Fig. 2 ). The radiative torque on the particle is given by the following expression: M r = p r (1 −0 ) ˆ t × S l ( ˆ t ·ˆ n )d S + 2 p r 0  S ˆ n ×l ( ˆ t ·ˆ n ) 2 d S, (6) where, as before, l is the vector from the particle centre of mass to the surface element d S . To test the validity of this geometric optics approximation, we compared the torque efficiency resulting from those expressions Figure 2. Fresnel reflection coefficients as a function of the incidence angle calculated for the assumed refractive index m = 1.6 + 0.2 i . The dashed lines correspond to the s and p polarizations of the incident radiation, and the solid line is the average of those two coefficients, which is the value used for 0 in equation (5). and those computed using the DDA code for oblate and prolate spheroids with the largest possible particle dimensions to make the geometric optics approximation valid. In all cases, a refractive index MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 Dynamics of cometary irregular particles 5145 Figure 3. X-component of the torque efficiency as a function of the angle of attack ξ(see Fig. 3 ) for oblate ( = 0.5 and = 0.125, left and middle hand) and prolate ( = 2.0, right hand) spheroidal particles of different ef fecti ve radii, as indicated. The refractive index is m = 1.6 + 0.2 i . The solid lines correspond to the geometric optics approximation described in Section 2.3, and the open circles joined by dashed lines correspond to DDA calculations. of m = 1.6 + 0.2 i is used. This value is appropriate for cometary particles composed by a mixture of silicates and a strongly absorbing component of carbonaceous and organic materials, and is in the range of values assumed by e.g. Markkanen et al. ( 2018 ) and Moreno et al. ( 2018 ) in their interpretation of the Rosetta/OSIRIS phase function measurements of the 67P coma particles. The torque efficiency vector is given by (Draine & Flatau 1994 ) as: Q M = kM πr 2 eff u rad , (7) where k = 2 π/ λ, and u rad = ¯ k r 2 h is the energy density due to solar irradiance at the heliocentric distance r h . The DDA calculations involve the use of dipole arrays obeying the constraint | m | kd ≤ 1, where d is the interdipole distance. This, in practice, limits the particle size to r eff ࣠ 5 μm for an incident wavelength of λ= 0.6 μm, owing to memory and CPU available resources. In our tests, we used oblate and prolate spheroids having axes ratios of = 0.5 and 2, with r eff = 5 μm, as well as a more extreme case of an oblate spheroid with = 0.125 and r eff = 7 μm. The number of dipoles needed to attain | m | kd ≈1 varied from 4 ×10 6 to 1.4 ×10 7 . The calculation involved between 8 and 20 h of CPU time per each particle orientation on a Dell R workstation with an Intel R Xeon R E5-1620 3.50GHz processor and 16 GB memory. In Fig. 3 , we display the DDA results together with our geometric optics approximation. The geometry of the problem is depicted in Fig. 3 , where the direction of the incident flux is along + Z. The X-component of the torque is depicted as a function of the angle of attack ξ. As it can be seen, the agreement between both calculations is good, even in the case of the oblate spheroid with large aspect ratio. This allows us to confidently use the geometric optics approximation as a fast algorithm to compute the radiative torques on the assumed absorbing and large particles in comparison with the incident wavelength. It is expected that for still larger particles that those of the comparison tests, the agreement would be even better, as the geometric optics regime would be fully attained. 2.4 Dynamical equations We used two reference frames, one centred on the particle centre of mass, denoted as ( xyz ), and another centred on the nucleus centre of mass, denoted by capital letters as ( XYZ ). The orientation of the particle with respect to frame ( XYZ ) is given by the classical Euler angles, ( φ, θ, ψ). Rotations of the particle are defined according to the (3,1,3) (or zxz ) Euler rotation sequence, also known as the x-convention (see e.g. Diebel 2006 ). The translational motion of the particle in the nucleus-centred reference frame is go v erned by the equation: m d d 2 r dt 2 = F a + F r + F gn + F gs , (8) where we include in the equation the nucleus gravity ( F gn ) and the solar gravity ( F gs ). The rotation of the particle is described in the ( xyz ) reference frame through the Euler equations. The particle axes are set to be coincident with its principal axes of inertia, so that the inertia tensor is diagonal and the Euler equations become: I xx dω x dt + ( I zz −I yy ) ω y ω z = M x (9) I yy dω y dt + ( I xx −I zz ) ω z ω x = M y (10) I zz dω z dt + ( I yy −I xx ) ω x ω y = M z , (11) where M = ( M x , M y , M z ) = M a + M r , ω x , ω y , ω z are the components of the angular velocity, and I xx , I yy , I zz , are the principal moments of inertia of the particle. To describe the linear and rotational motion of the particle, equations (8–11) in combination with the Euler kinematic equations must be numerically integrated. Ho we ver, in addition to numerical instabilities, there are difficulties when directly solving the Euler equations, such as the well-known Gimbal lock problem (see e.g. Zhao & van Wachem 2013 ). Instead, a numerical scheme based on the use of unit quaternions is much more accurate, and free of Gimbal lock. To this end, we followed the procedures given by Zhao & van Wachem ( 2013 ), which we detailed in the following subsection for completeness. 2.4.1 Quaternions: Algebra and relation to rotation matrices Quaternions were first introduced by Sir Hamilton (1844), and have been widely applied to solve mechanical problems since the XIX century. A quaternion q consists of a scalar part and a vector part, q = [ q o , q ] , and is defined as: q = q o + q 1 i + q 2 j + q 3 k , (12) where q o , q 1 , q 2 , q 3 are real numbers and ( i , j , k ) are unit vectors in the direction of the ( x , y , z) axes, respectively. The conjugate of a quaternion is given by q  = q o −q 1 i −q 2 j −q 3 k , while the norm is given by  q =  q 2 o + q 2 1 + q 2 2 + q 2 3 . A unit quaternion has norm of unity. The inverse of a quaternion is given by q −1 = q  /  q  , so that for a unit quaternion q one has q −1 = q  . The multiplication of two quaternions is defined as: t = pq = [ p o q o −pq , p o q + q o p + p ×q ] . (13) In rotation dynamics, a vector s can be rotated by applying a rotation matrix, which is a 3 ×3 matrix R where R T = R −1 and det( R ) = ±1. The rotated vector s  , is written by s  = R s . A coordinate (also called elemental) rotation is a rotation about a given MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 5146 F. Moreno et al. axis. The following three coordinate rotation matrices rotate a vector by an angle αabout the axis x , y , z: R x ( α) = ⎡ ⎣ 1 0 0 0 cos α−sin α 0 sin αcos α⎤ ⎦ R y ( α) = ⎡ ⎣ cos α0 sin α 0 1 0 −sin α0 cos α⎤ ⎦ R z ( α) = ⎡ ⎣ cos α−sin α0 sin αcos α0 0 0 1 ⎤ ⎦ . Any rotation of a vector s can be described by a sequence of three coordinate rotations. The most common rotation sequence is the (3,1,3), also known as the x-convention, that is the one that we have used in the modelling. The rotation matrix for a set of Euler angles ( φ, θ, ψ) in this (3,1,3) sequence is written as R zxz ( φ, θ, ψ) = R z ( φ) R x ( θ) R z ( ψ). Equi v alently, it can be shown that a vector s can be rotated by a unit quaternion q , resulting in a vector s  , by the expression (see e.g. Diebel 2006 ): s  = q s q −1 , (14) where the multiplication of a vector by a quaternion is simply the product of two quaternions where the vector is extended to a quaternion having zero scalar part and the same ijk components as the vector’s xyz components. Equation (14) implies a relationship between rotation matrices and unit quaternions. In the (3,1,3) sequence, it can be shown that a set of Euler angles ( φ, θ, ψ) can be converted to a quaternion using the equations (Diebel 2006 ): q o = cos ( φ/ 2) cos ( θ/ 2) cos ( ψ/ 2) −sin ( φ/ 2) cos ( θ/ 2) sin ( ψ/ 2) q 1 = cos ( φ/ 2) sin ( θ/ 2) cos ( ψ/ 2) + sin ( φ/ 2) sin ( θ/ 2) sin ( ψ/ 2) q 2 = cos ( φ/ 2) sin ( θ/ 2) sin ( ψ/ 2) −sin ( φ/ 2) sin ( θ/ 2) cos ( ψ/ 2) q 3 = cos ( φ/ 2) cos ( θ/ 2) sin ( ψ/ 2) + sin ( φ/ 2) cos ( θ/ 2) cos ( ψ/ 2) . The inverse transformation, quaternion to Euler angles, is given by the equations (Diebel 2006 ): φ= arctan 2 q 1 q 3 −2 q o q 2 2 q 2 q 3 + 2 q o q 1 θ= acos q 2 3 −q 2 2 −q 2 1 + q 2 o  ψ = arctan 2 q 1 q 3 + 2 q o q 2 −2 q 2 q 3 + 2 q o q 1 . 2.4.2 Numerical integration of the dynamical equations As stated before, to perform the numerical integration of the dynamical equations, we follow the predictor-corrector method posed by Zhao & van Wachem ( 2013 ). The particle coordinate axes are set first coincident with the comet (world) reference frame. Then, for the irregular particles, at time t = 0, we perform an initial random rotation of the particle through angles ( φo , θo , ψ o ). To perform the v alidation tests described belo w, we used oblate and prolate spheroids whose initial orientation is simply described by the angle of attack ξ(see Fig. 4 ). The particle is assumed at rest on the comet surface, and having zero angular velocity in the comet reference frame coordinates, i.e. V X = V Y = V Z = 0 and ω X = ω Y = ω Z = 0. F or conv enience, we assume that the initial coordinates of the particle are (0, 0, R n ). We assume that Figure 4. Initial position of a spheroidal particle in the comet reference frame. The direction of the gas and solar fluxes are along z -axis. The x -axis points towards the observer (right-handed system). ξis the angle of attack (see the text). the sun is placed at a distance Z = r h from the comet nucleus and located in the plane YZ . The particle is accelerated outwards by the gas drag in the + Z direction. In what follows, we will denote the vector components referred to the particle coordinate system with a superscript b (e.g. ω b or M b ) while no superscript is used for the variables expressed in the comet, or world reference frame. When indicating the current time level, we used a n subscript. The transformation of the torque and the angular velocity from the comet reference frame to the particle reference frame at a given time-step n is given by: ω b n = q −1 n ω n q n (15) M b n = q −1 n M n q n . (16) The Euler equations is then used to compute ˙ω b n (equations 8–10). Then, the angular velocities in the particle reference frame at time n + 1 4 and n + 1 2 are computed as: ω b n + 1 4 = ω b n + 1 4 ˙ω b n t (17) ω b n + 1 2 = ω b n + 1 2 ˙ω b n t. (18) The angular velocity in the comet reference frame at time n + 1 4 , ω n + 1 4 , is approximated by the quaternion q n as: ω n + 1 4 = q n ω b n + 1 4 q −1 n . (19) A prediction of the unit quaternion at time-step n + 1 2 , q  n + 1 2 is provided by the following equation: q  n + 1 2 =  cos  ω n + 1 4  t 4 , sin  ω n + 1 4  t 4 ω n + 1 4  ω n + 1 4   q n . (20) Then, the angular velocity in the comet reference frame at time n + 1 2 is obtained as: ω n + 1 2 = q  n + 1 2 ω b n + 1 2 q −1 n + 1 2 . (21) The torque at time n + 1 2 is given by: M n + 1 2 = q  n + 1 2 M b n + 1 2 q −1 n + 1 2 . (22) MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 Dynamics of cometary irregular particles 5147 Table 1. Model parameters used for validation with Ivanovski et al. ( 2017b ) calculations. Nucleus radius, R N [m] 2000 Nucleus mass, M N [kg] 10 13 Nucleus surface temperature, T N [K] 200 Gas composition H 2 O Gas specific heat ratio, γ1.33 Particle temperature, T d , [K] 200 The angular acceleration in the particle reference frame at time n + 1 2 , ˙ω b n + 1 2 can be then obtained from the Euler equation. The corrected quaternion at the next time-step n + 1 is given by: q n + 1 =  cos  ω n + 1 2  t 2 , sin  ω n + 1 2  t 2 ω n + 1 2  ω n + 1 2   q n . (23) And the angular velocity in the particle reference frame and in the comet reference frame at time-step n + 1 can finally be obtained as: ω b n + 1 = ω b n + ˙ω b n + 1 2 t (24) ω n + 1 = q n + 1 ω b n + 1 q −1 n + 1 . (25) At each time-step, we simultaneously compute the particle attitude, and the position and velocity of the particle with respect to the comet reference frame by equation (8), using the Euler method. The computational time-step must be much smaller than the inverse of the rotational frequency of the particle, but much larger than the gas-grain collisional time-scale. We used a time-dependent time-step given by t = 10 −3 2 π | ω| s. As a check to the validity of the computations, at each time-step we perform a correlation in the three spatial dimensions between the leftand right hand sides of the three Euler equations (9–11). We stop the e x ecution when the correlation coefficient becomes smaller than 0.99, indicating the presence of instabilities in the solution. This al w ays occur when the rotational frequencies are too high so that  t must be set to values even smaller than those given by the abo v e e xpression, making the problem intractable in practice. This problem mostly appears when dealing with small ( r eff ࣠ 10 μm) particle sizes as shown later in Section 3. 2.4.3 Sample execution of the code: Validation and study of the effect of the combined gas plus radiation torques We have validated our code via comparisons with results obtained by Ivanovski et al. ( 2017b ), which refer to gravitational attraction of the nucleus, and aerodynamic force and torque on spheroidal particles. For these comparisons, the model parameters are those of Table 1 by Ivanovski et al. ( 2017b ), that we reproduce here for completeness in Table 1 . The spheroidal particles experience an oscillatory behaviour first, and eventually, at time t rot , the particle starts to display full rotation. The asymptotic rotation frequency is denoted by ν∞ . The particle is accelerated outwards from the comet surface, reaching eventually a terminal velocity, V ∞ . The first test was made with an oblate spheroid of axial ratio a / b = 0.5, ef fecti ve radius of 1 mm, and density of ρ= 100 kg m −3 , placed on the comet surface at an angle of attack of 45 ◦, subjected to a water production rate of Q g = 10 28 molecules s −1 . This corresponds to case #b03g3d1ob of table B.9 by Ivanosvki et al. ( 2017b ). Our calculations yielded V ∞ = 15.2 m s −1 , an asymptotic frequency of ν∞ = 0.107 s −1 , and a t rot = 535 s, in excellent agreement with Ivanovski et al’s results ( V ∞ = 15 m s −1 , ν∞ = 0.1, and t rot = 534 s). Fig. 5 displays the Figure 5. Rotational frequency and velocity as a function of time for an oblate spheroid having = 0.5 and ef fecti ve radius r eff = 1 mm. The physical parameters for the comet and gas environment are given in Table 1 . Figure 6. The modules of the aerodynamic (black line) and radiative (red line) torques as a function of time and distance to the nucleus for the oblate spheroid of = 0.5 and r eff = 1 mm described in the text. The physical parameters for the comet and gas environment are given in Table 1 . rotational frequency and velocity as a function of time for such particle. Additional tests were made for a variety of conditions, including different production rates, densities, angles of attack, and shape models (prolate versus oblate). In all cases, the agreement with the computations of Ivanovski et al. ( 2017b ) were excellent. The next step is to show some results aimed at illustrating the rele v ance of including the combined effect of the gas drag and radiative forces and torques. Fig. 6 shows the modules of the aerodynamic and radiative torques on the oblate spheroidal particle described abo v e. In this graph, we see that the aerodynamic torques clearly dominate until t ∼7000 s or to a distance of the nucleus centre of R ∼100 km, where the radiative torques start to compete with the aerodynamic torques. But even at larger distances, the effect of the radiative torques on such particle remains negligible. Fig. 7 displays the behaviour of the frequency and particle velocity for an integration time of 2.4 ×10 4 s (i.e. ∼360 km from the nucleus), where we see that the rotation frequency remains unaltered. Only a very slightly decrease of the velocity, as a consequence of the radiative force, which is opposite to the gas flux, is noticed. For oblate spheroids of smaller ef fecti ve radius of r eff = 10 μm, the rotation frequencies encountered are higher, as expected, but the effect of radiation torques MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 5148 F. Moreno et al. Figure 7. Long-term evolution of the rotational frequency and velocity for an oblate spheroid having = 0.5 and ef fecti ve radius r eff = 1 mm. The physical parameters for the comet and gas environment are given in Table 1 . Figure 8. Long-term evolution of the rotational frequency and velocity for an oblate spheroid having = 0.5 and ef fecti ve radius r eff = 10 μm. The physical parameters for the comet and gas environment are given in Table 1 . remain negligible, as in the larger particle size regime. Fig. 8 shows the rotation frequency and velocity for such small particle, were we see that the asymptotic frequency is ∼8 Hz, and the velocity, after reaching a maximum value near 140 m s −1 decrease because of radiation pressure to ∼130 m s −1 at the latest position recorded ( ∼500 km from the nucleus at t = 3690 s). We underline that the effect of the radiative torque on the particles, although being negligible compared to the expanding gas torques in the nucleus vicinity, they might become important on longer timescales. Ho we ver, on long time-scales other processes, such as the Yarko vsk y-O’K eefe-Radzie vskii-Paddack ef fect (see e.g. Rubincam 2000 ), or the Poynting–Robertson drag, become important as well. The building-up of a model taking into account these other forces is well beyond the scope of this work. 3 RESULTS OF THE DYNAMICAL MODELLING The purpose of this section is first to show the results of the dynamical models applied to the three irregularly shaped particles of Fig. 1 under different gas and solar radiation environments, and for a variety of effective sizes. Then, in a dedicated subsection, we build-up a larger data base of particles in a wide axial ratio distribution to study their dynamical evolution in order to asses more firmly the dynamical parameters. In order to keep the number of free parameters of the model to a minimum, we will restrict ourselves to the results rele v ant to the physical environment of Table 2. Model parameters appropriate for comet 67P/Churyumov– Gerasimenko. Parameter Value Reference Nucleus density [kg m −3 ] 533 P ¨ atzold et al. ( 2016 ) Nucleus volume [km 3 ] 18.7 P ¨ atzold et al. ( 2016 ) Nucleus eff. radius [m] 1650 Q g [#s −1 ] ( r h = 1.24 au) 3 ×10 28 Hansen et al. ( 2016 ) Q g [#s −1 ] ( r h = 2.00 au) 10 27 Hansen et al. ( 2016 ) Particle eff. radius [m] 10 −5 to 10 −3 Particle density [kg m −3 ] 800 Fulle et al. ( 2016 ) P article refractiv e inde x m = 1.6 + 0.2 i Table 3. Dynamical parameters of r eff = 1 mm particles at two heliocentric distances. Particle properties at perihelion ( r h = 1.24 au) Shape DT L fi n a l lat. νfinal V final Model (deg) (deg) (Hz) (m s −1 ) Begonia 12 ±4 –1 ±25 1.6 ±0.7 10.0 ±0.2 Glasenappia 23 ±12 + 4 ±32 2.7 ±1.1 10.0 ±0.3 Aurora 2.1 ±0.7 + 3 ±25 0.9 ±0.5 9.34 ±0.05 Particle properties at r h = 2 au Shape DT L fi n a l lat. νfinal V final Model (deg) (deg) (Hz) (m s −1 ) Begonia 13 ±3 –5 ±23 0.3 ±0.2 1.60 ±0.05 Glasenappia 24 ±10 + 5 ±32 0.5 ±0.3 1.60 ±0.05 Aurora 2.2 ±0.6 + 6 ±22 0.15 ±0.07 1.45 ±0.01 comet 67P/Churyumov–Gerasimenko, as it has been subject of much interest because of the successful European Space Agency mission Rosetta . The comet environmental parameters are those show in Table 2 . The remaining parameters (dust and surface temperatures, T d , T S , and γ) are assumed as in Table 1 . In order to establish the ensemble properties for each particle shape when subjected to a certain environment, we place the selected particle shape at the nucleus surface and then perform a random rotation sequence in the Euler angles φ, ψ, and θaccording to φ= 2 πr 1 , ψ = 2 πr 2 , and θ= acos (1 −2 r 3 ), where r 1 , r 2 , and r 3 are random numbers in the [0,1] interval. We then followed the trajectory of each particle and record its attitude up to a distance from the nucleus of 50 km, where the gas drag becomes negligible. The procedure is then repeated for 100 initial positions of each particle shape. Ho we ver, although our initial purpose was to accomplish this task for particles having ef fecti ve radii of 10 μm, 100 μm, 1 mm, and 1 cm, the procedure could not be fully completed for the two smallest sizes (10 and 100 μm) owing to the presence of instabilities in the dynamical evolution that appear for most initial conditions, as explained below. For each particle shape and ef fecti ve radius (1 mm and 1 cm), and for the two heliocentric distances selected, we calculated the degree of tumbling, DT , defined as the mean angle between the angular momentum vector and the spin axis along the integration, and the final direction of the angular momentum vector L fi n a l , the final rotation frequency ( νfi n a l =  ν2 x + ν2 y + ν2 z , where ν= ω/2 π), and the final velocity V final at the ending time of integration, when the particles reach a distance of 50 km from the nucleus centre. The two first parameters, DT , and L fi n a l are computed mainly for the sake of comparison with the calculations by ˇ Capek ( 2014 ). Tables 3 and 4 give the obtained parameters for ef fecti ve radii of 1 mm and 1 cm. In the former case, the parameters pertain to two MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 Dynamics of cometary irregular particles 5149 Table 4. Dynamical parameters of r eff = 1 cm particles at perihelion ( r h = 1.24 au). Shape DT L fi n a l lat. νfinal V final Model (deg) (deg) (Hz) (m s −1 ) Begonia 12 ±4 0 ±21 0.09 ±0.06 3.06 ±0.07 Glasenappia 26 ±10 + 6 ±33 0.2 ±0.1 3.10 ±0.09 Aurora 2.2 ±0.5 + 1 ±30 0.04 ±0.02 2.85 ±0.02 heliocentric distances, while in the latter case, only the parameters at perihelion could be obtained, as the modelled gas drag is not enough to lift 1 cm particles at r h = 2 au. Concerning the degree of tumbling, we conclude that it depends on the shape model used, but neither on the size nor on the heliocentric distance. The rounded shape model, Aurora , shows the smallest DT ∼2 ◦, while the elongated shape show the largest DT ∼25 ◦. Interestingly, the flattened shape, Begonia , shows the same value DT ∼12 ◦as that found by ˇ Capek ( 2014 ). Considering all three shape models together, the median of results for r eff = 1 mm at perihelion gives DT = 12 ±9 ◦, while a very similar value for r eff = 1 cm ( DT = 13 ±9 ◦) is found. At r h = 2 au, the median of the three shape models for r eff = 1 mm gives again a similar result ( DT = 13 ±8 ◦). Another interesting result is the fact that the latitude of the angular momentum vectors at the end of the integration point to a latitude interval al w ays centred around 0 ◦, i.e. perpendicular to the gas flux direction. This is also consistent with ˇ Capek ( 2014 ) findings. The final frequencies encountered are clearly dependent on the particle size and shape, and on the heliocentric distance. The largest frequencies are found for the most elongated shape, Glasenappia , which vary between 0.2 Hz for 1 cm particles and 2.7 Hz for 1 mm particles at perihelion. For the most rounded shape model, Aurora , notably smaller frequencies are found, moving in the 0.04–0.9 Hz range. The flattened particle rotates at intermediate frequencies between the other two model shapes. Overall, the values found for the frequencies of mm to cm sized particles are very consistent with the those values found by Frattin et al. ( 2021 ) in their analysis of a large number of single particle tracks from OSIRIS camera aboard Rosetta . From the light curves, Frattin et al. ( 2021 ) derived frequencies in the 0.24– 3.6 Hz range, with a strong maximum near the lower limit. Ho we ver, o wing to the limitations of the measurements, the most probable frequency is likely below the sampled lower limit, which allows us to constrain the particle size to be most probably in the cm or higher range. In addition, many of the particle tracks seen by OSIRIS have associated light curves corresponding to complex rotational motion, as found here from the irregularly shaped particles dynamics. Owing to the irregularity of the particles, the rotation frequencies are al w ays considerably higher than symmetric particles of the same ef fecti ve radii and similar aspect ratios. In Fig. 9 , a fe w typical examples of the evolution of the frequency components and the rotational energy of particles of different sizes and shapes at the conditions of 67P perihelion are shown. In contrast with those symmetric particles, the rotation is chaotic right after ejection from the surface, showing a complex rotational motion. After some distance from the nucleus, the frequency components evolve e ventually to wards a more regular pattern, and the rotational energies tend to asymptotic values (see Fig. 9 ). A similar behaviour was found by ˇ Capek ( 2014 ) in his analysis of rotating meteoroids. ˇ Capek ( 2014 ) also found a relationship between median spin frequencies and velocities for the cases of diffuse and specular reflection of gas molecules on the particle surface. For diffuse reflection, he found ¯ν≈2 ×10 −3 V fi n a l D −0 . 88 , where D is the particle diameter. For particles of 1 cm radius and V final ∼3 m s −1 , ¯ν∼0.19 Hz, which is of the same order of the frequencies shown in Table 4 , particularly for the Glasenappia shape model particle, and also in good agreement with the frequencies obtained for flattened and elongated particles having wide axial ratio distribution as described below in Section 3.1 (see Table 6 ). For smaller (1 mm radius) particles, the predicted frequencies at perihelion and r h = 2 au are ∼4.5 and ∼0.75 Hz, respectively, which compare well, although are bit higher than those shown in Table 3 . As noted abo v e, when the particle ef fecti ve radius is set down to 100 μm, in quite a few cases there appear numerical instabilities that prevent us to perform any meaningful statistics. The situation becomes worse for still smaller particles. Successful examples of the dynamical evolution of a r eff = 100 μm Begonia shape model particle and of a r ef f = 10 μm Glasenappia shape model particle are given in Fig. 9 (lower panels). For particles having r eff = 100 μm, rotation frequencies of order 20–40 Hz can be reached, while for ∼10 μm particles, frequencies as high as 500–1000 Hz are observed. In both cases, as for the larger particles, after a period of time of chaotic rotation close to surface, the frequencies evolve to a more smooth behaviour as the gas drag becomes less and less important. The high frequencies attained imply that the time-dependent time-step  t had to be as small as ∼10 −7 s in order to keep the dynamical solution stable during the integration. But even so, numerical instabilities are seen to appear for most initial attitudes of the particles at a few km from the surface, suggesting that the integration time-step should be still smaller, so that the problem becomes intractable in practice. The final velocities, V final are found to depend on heliocentric distance and size, as expected, but only very slightly on the particle shape. The maximum values are found for the elongated and flattened shapes, with values very close between them, and higher than those shown by the rounded shape in the 7–9 per cent interval only. This behaviour is similar to that shown by spheroidal particles when compared with spheres of same ef fecti ve radius (Iv anovski et al. 2017a ). As stated abo v e, the numerical simulation for fast-rotating particles requires long computational time. On the other side, if we are interested in an assessment of order of magnitude, it is possible to scale the available numerical data to other sizes instead of time-consuming numerical simulations. For the scaling, we use the relation from Ivanovski et al. ( 2021 ) which is based on the approach proposed in Zakharov et al. ( 2018 ) and further developed in Zakharov et al. ( 2021 ). This approach uses a set of universal, dimensionless parameters, which characterize the dust motion in the inner cometary coma and allows one to reveal dust flows similarities. Table 5 shows comparison of scaled data with the available numerical results for the smallest sizes analysed, 10 and 100 μm. The scaled velocities are in a very good agreement with the ones calculated numerically. Although the scaled frequencies differ from the numerical values up to 3–6 times, they still provide a reasonable order-of-magnitude estimate, taking into account the uncertainties in the computed frequencies. The scaled frequencies for r eff down to 1 μm are also given in Table 5 . This reveals that rotation frequencies in excess of 1000 Hz can be attained for such small particles. A high rotation frequency could led to particle disruption if centrifugal forces exceed the tensile strength. Ho we ver, follo wing estimates of the critical rotation period by Davidsson ( 1999 ), the tensile strength of the material composing the 1 μm particle rotating at 1000 Hz MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022 5150 F. Moreno et al. Figure 9. Examples of the evolution of the rotational frequency components νx , νy , and νz (red, green, and blue lines, respectively, left hand side y -axes), and the corresponding rotational energy (black lines, right hand side y -axes), as a function of distance to the nucleus centre, for various model particles with different sizes as indicated. The left hand correspond to a Begonia -shaped particle, while the right hand pertain to a Glasenappia -shaped particle. The environment conditions are those of comet 67P at perihelion (see Table 2 ). Table 5. Comparison of scaled and computed values of terminal velocity V final and frequency νfinal . Shape V final (m s −1 ) νfinal [Hz] Scaled Computed Scaled Computed Scaling from r eff = 1 mm (Tab.3) to r eff = 100 μm Begonia 31.62 31.50 15.62 97.50 Scaling from r eff = 1 mm (Tab.3) to r eff = 10 μm Glasenappia 100.00 98.20 243.57 780.00 Scaling from r eff = 1 mm (Tab.3) to r eff = 1 μm Begonia 316.23 –1067.13 – Glasenappia 316.23 –1800.78 – Aurora 295.36 –620.24 – should be ࣠ 10 −3 Pa to breakup, which is several orders of magnitude below the current estimates of tensile strengths for cometary solid particles (G ¨ uttler et al. 2019 ). Another aspect of the dynamics of those rotating particles which is important to address is whether their associated scattering phase function would show a strong backscattering enhancement with a minimum near 90 ◦phase angle, based on a purely geometric effect. This minimum has been clearly found by Bertini et al. ( 2017 ) in 67P dust from Rosetta /OSIRIS images, and the backscattering effect has also been seen in several other comets (see the compilation of observations in Bertini et al. 2017 ). This phase function shape has been interpreted with a complex light scattering model considering the presence of large particles constituted by densely packed submicrometer-sized grains of organic material and larger micrometer-sized grains of silicate composition (Markkanen et al. 2018 ). In a previous work (Moreno et al. 2018 ), we provided an alternative interpretation based on the fact that a combination of large prolate and oblate shaped particles oriented in such a way that their shorter axes point toward the sun would easily explain the shape of the observed phase function by OSIRIS as well, without being strongly dependent on the composition or the structure of the particles. Laboratory scattering measurements performed by Mu ˜ noz et al. ( 2020 ) with large oblate-shaped porous and absorbing particles confirm the theoretical predictions. With our dynamical model, we can perform the calculations that allows us to check whether this MNRAS 510, 5142–5153 (2022) Downloaded from https://academic.oup.com/mnras/article/510/4/5142/6484808 by Inst. Astrofisica Andalucia CSIC user on 04 May 2022