Issue 
A&A
Volume 643, November 2020



Article Number  A16  
Number of page(s)  7  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202039092  
Published online  27 October 2020 
Thermophysical model for icy cometary dust particles
^{1}
Max Planck Institute for Solar System Research,
JustusvonLiebigWeg 3,
37077
Göttingen, Germany
email: markkanen@mps.mpg.de
^{2}
Institut für Geophysik und extraterrestrische Physik, Technische Universität Braunschweig,
Mendelssohnstr. 3,
38106
Braunschweig, Germany
Received:
3
August
2020
Accepted:
10
September
2020
Context. Cometary dust particles are subjected to various forces after being lifted off the nucleus. These forces define the dynamics of dust, trajectories, alignment, and fragmentation, which, in turn, have a significant effect on the particle distribution in the coma.
Aims. We develop a numerical thermophysical model that is applicable to icy cometary dust to study the forces attributed to the sublimation of ice.
Methods. We extended the recently introduced synoptic model for icefree dust particles to icecontaining dust. We introduced an additional source term to the energy balance equation accounting for the heat of sublimation and condensation. We use the direct simulation Monte Carlo approach with the dusty gas model to solve the mass balance equation and the energy balance equation simultaneously.
Results. The numerical tests show that the proposed method can be applied for dust particles covering the size range from tens of microns to centimetres with a moderate computational cost. We predict that for an assumed ice volume fraction of 0.05, particles with a radius, r ≫ 1 mm, at 1.35 AU, may disintegrate into mmsized fragments due to internal pressure buildup. Particles with r < 1 cm lose their ice content within minutes. Hence, we expect that only particles with r > 1 cm may demonstrate sustained sublimation and the resulting outgassing forces.
Key words: radiative transfer / diffusion / methods: numerical / comets: general
© J. Markkanen and J. Agarwal 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1 Introduction
As comets approach the Sun, they start to eject gas and dust from their nuclei. The ejection of dust is mainly driven by the sublimation of volatile ices. Thermal modelling suggests that the sublimation of volatiles such as CO and CO_{2} ices can lift large particles containing water ice off the surface, whereas sublimation of water ice can only lift smaller, icefree particles (Gundlach et al. 2020). In addition, other mechanisms, such as the crystallization of amorphous ices and impacts, may be capable of ejecting smaller icy particles. The icecontaining particles ejected from the surface may behave like minicomets whose trajectories are affected by the recoil force due to anisotropic outgassing of water vapour. Indeed, particles whose trajectories cannot solely be explained by gravity, gas drag, and radiation pressure have been detected for the comets 103P/Hartley 2 (Kelley et al. 2013, 2015) and 67P/ChuryumovGerasimenko (Agarwal et al. 2016). Furthermore, the sublimation of ices can build up pressure in the interior of a particle, which may lead to their disintegration while they are travelling from the inner to outer coma.
Understanding the outgassing of dust particles requires rigorous multiphysical modelling. Typical thermophysical models applied to the cometary nuclei, however, introduce various simplifications which make them inapplicable to small dust particles. More specifically, the models are often one dimensional, they assume wavelengthindependent surface absorption and emission of electromagnetic radiation, the radiative heat transfer is treated as a onedimensional transport by defining an effective heat transfer coefficient, and the gas flow is treated as the Knudsen flow (Huebner et al. 2006; Prialnik 2004; Gundlach et al. 2020). For small particles, these approximations may introduce significant errors. The thermal evolution of small cometary dust particles is often modelled by assuming isothermal spherical particles and using the analytical LorenzMie solution for the absorbed and emitted radiation and assuming that the sublimation only occurs on the surface of the particle (Lichtenegger & Kömle 1991; Gicquel et al. 2012, 2016). This is a valid assumption for small compact particles, but inside porous particles larger than tens of micrometres, temperature and pressure gradients may occur in a typical coma environment. Thus, a more realistic numerical model is needed to bridge the gap between small and very large particles.
We recently introduced a numerical method for analysing the thermal properties of icefree particles by treating the radiative heat transfer with the radiative transfer with the reciprocal transactions framework (Markkanen & Agarwal 2019). Here, we extend the method for icecontaining particles by introducing an additional source term to the energy equation that accounts for the latent heat of sublimation and condensation. Because this term depends on the temperature and pressure, we solve the energy balance equation simultaneously with the mass balance equation. We apply the finiteelement method (FEM) for the energy balance equation and the direct simulation Monte Carlo (DSMC) method with the dusty gas model for the mass transport.
The article is organized as follows. In Sect. 2, we present the governing equations describing the physics of the problem. The numerical methods we used are explained in Sect. 3. Section 4 introduces the particle model. In Sect. 5, we present the validation of the developed method against the continuous gas flow model and apply the method to study thermal properties of icy cometary dust particles at 1.35 AU. Finally, we present our conclusions in Sect. 6.
2 Governing equations
Recently, we introduced a thermal model for the aggregated icefree particles composed of submicrometresized dust grains in Markkanen & Agarwal (2019). Here, we will extend the method for dust particles containing ice inclusions. We assume that the dust and icegrains are intimately mixed in the microscopic scale. We also assume local equilibrium, in the case of which the solid and gas phases have the same temperature locally and the energy exchange between the solid and gas phases is neglected.
2.1 Energy balance
The energy balance equation for the solid phase is written as (1)
where c_{p} is the specific heat capacity, ρ_{s} is the density, T is the temperature, and κ is the conductive heat transfer coefficient. The righthand side includes the volumetric source terms, the absorbed solar and thermally emitted and reabsorbed radiation Q_{r}, and the latent heat of sublimation and condensation of ices Q_{l}. The above values are macroscopic, that is, they are averaged over a small volume element larger than the microstructure of dust and ice grains. We also assume that the density and specific heat capacity of gas are many orders of magnitude lower than those of the solid phase and ignore the heat transport in the gas phase.
On the particle surface, we impose the Neumann boundary condition for the conductive flux as n ⋅ κ∇T = 0, where n is the outer unit normal vector. Thus, only the radiative heat flow can cross the boundary as the convective flow is neglected.
Energy related to the sublimation and condensation processes can be written as (2)
where q is the sublimation rate and H is the latent heat. By defining the saturation vapour pressure as (3)
and using the experimentally obtained coefficients, A and B (Fanale & Salvail 1984), the heat of sublimation is obtained via the ClausiusClapeyron relation and given by (4)
in which R_{g} is the universal gas constant and μ is the molar mass of the gas.
The production rate is obtained via the Hertz–Knudsen formula as (5)
where S is the surfacetovolume ratio in the microscale, and P is the pressure in the micropores.
2.2 Massbalance
The mass balance equation for the gas and solid phases is written as
where ρ_{g} and ρ_{s} are the gas and solid phase densities, respectively, and J is the gas flux density. The macroscopic thermodynamic properties of gas such as the density, velocity, pressure, and flux emerge from the stochastic microscopic description of the state of gas given by the Boltzmann equation: (8)
where f(t, r, v) is the time (t) position (r) and velocity (v) dependent distribution density function, f_{q} is the source function, and C(f, f) is the collision operator. Formally, the collision operator reads (9)
where σ is the kernel that describes scattering, the primed functions are the post collision distribution densities, and the subscript * denotes the other particle in the collision pair.
The macroscopic number density n can be computed by integrating the distribution density function over the velocity space as (10)
We apply the DSMC method to solve Eq. (8), as described in Sect. 3.
3 Numerical solution
We discretize the domain of interest Ω with tetrahedral elements. Then we employ the FEM for the spatial dimension of the energy balance Eq. (1) with the nodal testing w^{n} and basis functions u^{m}. For the temporal dimension we will use the central finite difference formula (12)
and interpolate the other terms at as (13)
This results into the Crank–Nickolson scheme for the unknown coefficient vector x_{t+1} given by (14)
where the mass and the stiffness matrices are defined as (15) (16)
Equation (14) is strongly nonlinear because of the force term. Hence, we use an iterative method with underrelaxation to solve the unknown coefficient vector, x_{t+1}.
We compute the absorbed solar and thermally emitted radiation, Q_{r}, by using the radiative transfer with reciprocal transactions framework, as described by Markkanen & Agarwal (2019). Energy related to sublimation and condensation requires evaluation of the gas production rate, q, which depends on the temperature, T, and pressure, P. Thus, we need to solve the mass balance equation simultaneously with the energy balance equation. We apply the DSMC method, and its solution thus satisfies the Boltzmann equation, to solve the gas flow by assuming that the flow is a free molecular flow, that is, the mean free path of moleculemolecule collisions is much larger than the mean free path of moleculegrain collisions. The grains are assumed to be much heavier than the gas molecules and, thus, they are stationary, giving us the socalled dusty gas model. We consider the moleculegrain collisions stochastically, as in Ahmadian et al. (2019).
To compute the source function f_{q} in the DSMC, we introduce new (or remove the existing) DSMC molecules for each tetrahedral element with the index, n, based on the production rate, q^{n}, at any given time step, Δt_{DSMC}, with weight, q^{n}Δt_{DSMC}. The velocity of the new DSMC molecule is drawn from the local solidphase temperaturedependent MaxwellBoltzmann distribution. Then the trajectories of the DSMC molecules are traced within a time step. To account for the collision operator, C, the moleculegrain scattering distance is drawn from the exponential distribution with the mean free path, l_{DSMC}, and a new velocity for the gas molecule is drawn the MaxwellBoltzmann distribution defined by the local temperature. Finally, the macroscopic thermodynamical properties, for instance, density, pressure, and flux are sampled by summing up the DSMC molecules in each tetrahedron and averaging them over the time frame. The process is repeated until the steady state or the energy equation time step is reached.
The characteristic time scales for energy and mass transport can be very different. This means that we need to use different time steps (τ, Δt_{DSMC}) for the energy and mass transport to get a sufficiently fast numerical method. The DSMC time step, Δt_{DSMC}, must be small enough such that the DSMC molecules cannot travel across a finite cell within a single time step, whereas τ is limited by the nonlinearity and the convergence of the iterative method used to solve the energy balance equation. In practice, for cometary dust applications, the DSMC time step must be many orders of magnitudes smaller than the energy balance equation time step. Often, the gas flow reaches the steady state within the energy balance equation time step. In such a case, we stop the DSMC simulation and extrapolate the total gas production rate and the remaining ice mass at t + τ. The extrapolation does not conserve ice and gas mass exactly when a tetrahedron runs out of ice within the time step τ as the ice mass cannot be negative. To avoid the problem, we run the DSMC again by setting the ice mass to zero for the tetrahedron that has run out of ice and calculate the production rates using the total mass loss at the given tetrahedron as (18)
where m_{t} is the ice mass at time t.
4 Particle model
We used a particle microstructure model derived by fitting the phase function of the coma of 67P/ChuryumovGerasimenko in Markkanen et al. (2018). The heat conduction coefficient was derived by fitting the superheating phase function of the same comet in Markkanen & Agarwal (2019).
The microstructure consists of submicrometresized organic monomers and micrometresized silicate monomers randomly deposited inside a particle with the porosity, Φ = 0.6. The effective heat capacity, c_{p} = 750 J kg^{−1} K^{−1}, and the effective density, ρ_{s} = 1000 kg m^{−3}, are estimated from the porosity Φ and the silicatetoorganic volumetric ratio of 1/5. The size distributions of monomers follow a differential power law with the index of − 3, and the minimum and maximum cutoff limits are a_{min} = 65 nm and a_{max} = 125 nm for the organic and a_{min} = 0.6 μm and a_{max} = 1.3 μm for the silicate monomers. The refractive index for the organic monomers is approximated to be the same as amorphous carbon and is taken from Jäger et al. (1998). For the silicate monomers, the refractive index is taken from Dorschner et al. (1995). Here, we assume that ice is uniformly deposited in the micropores and the amount of ice is so low that it does not have an effect on the scattering properties. The microstructure model gives rise to the incoherent scattering properties of the volume elements which are used to evaluate the radiative part of the energy equation as described by Markkanen & Agarwal (2019). In addition, the gas mean free path, l_{DSMC}, and the surface to volume ratio, S, are derived from the microstructure model. The parameters used in the model are presented in Table 1.
The macrostructure is described by the tetrahedral mesh in which each tetrahedron can be assigned to different microstructure, effective material properties c_{p}, ρ_{s}, κ, and ice content. To generate the macrostructure model, we applied the hierarchical Voronoi partitioning algorithm with the parameters, N_{1} = 100, N_{2} = 1000, , and from Markkanen et al. (2015). This creates particles with the macroporosity of 0.6 and the mean unit Voronoi cell size of one tenth of the particle size. Figure 1 shows an example computational mesh generated by the twolevel hierarchical Voronoi algorithm. In the remainder of the paper, the size of the particle refers to the radius, r, of the smallest sphere containing the particle if not defined otherwise.
5 Results
5.1 Validation
To validate the DSMC implementation, we compared the results with the continuous gas flow model in porous media. The Fick’s law states that the gas flux J is proportional to the gas density ρ_{g} gradient as (19)
where D is the diffusion coefficient. Here, the gas density ρ_{g} is defined asthe gas mass per total volume including dust grains. Using the continuity equation, we have (20)
Writing the gas density in terms of pressure inside the micropores and using the ideal gas assumption, (21)
leads to the following equation for the gas density (22)
The diffusion coefficient for the random walk in three dimensions with the Maxwell–Boltzmann mean velocity v_{th} is given by (23)
Furthermore, by assuming the half MaxwellBoltzmann distribution above the particle interface, giving the normal mean velocity of v_{n} = v_{th}∕4 (Huebner & Markiewicz 2000), the Neumann boundary condition is given as (24)
Equation (22), supplemented with the boundary condition (24) can be solved with the FEM analogously to the FEM solution of the energy balance Eq. (14) by changing the coefficients and the righthand side accordingly and introducing the boundary integral for the Neumann boundary condition (which is zero for the energy equation).
First, we considered the simplest possible case; an isothermal T = 200 K spherical particle with the radius r = 0.1 mm and S = 10 000 m^{−1} in a steady state. No other energy sources were included. We solved the problem using the DSMC and the FEM. Two different DSMC solutions were computed, namely, DSMC 1 and DSMC 2, the latter having twice the number of the DSMC molecules than the former. Figure 2 shows the pressure sampled at the barycentre of each tetrahedron as a function of the distance from the centre of the spherical particle. The solutions are in good agreement with each other but the DSMC solutions have more noise than the FEM solution. The DSMC noise decreases with the increasing number of the DSMC molecules, as expected.
It is interesting to note that when S is constant, Eq. (22) is reduced to the Helmholtz equation in a steady state with the wavenumber squared given by (25)
The fundamental solution, that is, the Green’s function for the Helmholtz equation is written as (26)
where r and r′ are the source and observation points, respectively, and i is the imaginary unit. Thus, with the imaginary k the solution is an exponentially decaying evanescent wave, and the pressure drop near the surface can be very steep when SD^{−1} is large. This also means that in order to solve Eq. (22) with the standard continuous Galerkin FEM, an extremely fine mesh is needed at the sublimation front, which can make computational time prohibitively long.
Next, we simulated an irregularlyshaped dust particle under solar radiation at 1.35 AU. The initial temperature of the particle was 160 K and the volumetric ice content was 0.05. The surfacetovolume ratio parameter was set to S = 1 m^{−1} to allow for an efficient FEM solution. The volume equivalent radius of the particle was 0.5 mm and it was spinning 1 rpm with the spin axis pointing perpendicular to the direction of the Sun. The solar illumination was switched on at t = 0 s. The time evolution of the average temperature, gas production rate, maximum pressure averaged over a tetrahedron inside the particle and the ice volume fraction are plotted in Fig. 3 computed by using the DSMC and FEM for the gas transport. We observed an excellent agreement between the two different methods. We note, however, that the boundary condition used in the FEM does not account for the gas flow back to the particle, whereas it is taken into account in the DSMC method. The gas production rate and maximum pressure peak at t ≈ 200 s and then they decrease. This happens because sublimation creates an insulating dry dust layer on the particle’s surface that dampens the energy transport into the interior of the particle.
Simulation parameters.
Fig. 1 Crosscut of an example computational mesh generated by the twolevel Voronoi partitioning. The yellow tetrahedra correspond to the particle with a specified microstructure and the blue ones correspond to free space. 
Fig. 2 Pressure as a function of distance from the centre of the sphere calculated by the DSMC and FEM. The DSMC 2 solution has twice the number of the DSMC molecules than the DSMC 1 solution to demonstrate convergence. 
Fig. 3 Comparisons between thermal modelling results obtained by using the DSMC and FEM solutions for gas transport. The test particle has an irregular shape, with the volume equivalent radius of 0.5 mm. The surfacetovolume ratio was set to S = 1 m^{−1}. 
5.2 Application to cometary dust
Next, we studied thermal properties of dust in a cometary coma at 1.35 AU. We assumed that the ejected dust particles have temperature T_{0} = 160 K and contain5 vol. of water ice uniformly distributed into the micropores at t = 0. Thus, the dust particle was assumed to originate below the hot surface layer in order to contain ice. Here, we do not account for the ejection mechanism as our goal is to study the thermal evolution of icy dust particles once they have been ejected from the nucleus and exposed to the direct sunlight. In the simulation, we assumed that the particles are rotating around the axis perpendicular to the solar direction with the angular speed of 1 rpm. The results were not averaged over an ensemble of particles due to computational time restrictions.
Figure 4 plots the ice volume fraction as a function of time for 0.1 mm, 1 mm, and 10 mmsized particles. We observed that the particleswith r = 0.1 mm run out of ice in approximately a few tens of seconds after the ejection, r = 1 mm in ten minutes, and r = 10 mm particles can be extrapolated to stop subliming after a few hours. We expect that either they run out of ice or create an insulation surface layer, completely dumping sublimation. The time integration was stopped after 20 days of computing using 24 cores for the particle with r = 10 mm. The smaller particles took less than a week each with 24 CPUs to complete the simulation. The biggest bottleneck in the simulations of large particles is the nonlinearity of Eq. (14). A very small time step τ is required to solve the nonlinear equation prohibiting long time evolution simulations. Thus, using a more sophisticated nonlinear iterative solver or finding an efficient preconditioner would be an interesting topic for future research. Also, the DSMC simulation becomes slower with the increasing particle size as the diffusion time is proportional to r^{2}. Thus, for large particles, a continuous gas transport model would be preferred, assuming that the stability problem appearing for high S due to the exponentially decaying pressure profiles can be solved.
The sublimation of ices increases pressure inside the particles as demonstrated in Fig. 5. If the pressure is high enough, the particle may disintegrate. An experimentally verified model (Skorov & Blum 2012; Blum et al. 2014; Brisset et al. 2016) for the effective tensile strength of a particle made of aggregates is given by Pa where r_{a} is the radius of aggregates. In our particle model, these aggregates can be considered to correspond the Voronoi cells with r_{a} = r∕10. The tensile strengths for the particles made of aggregates with r_{a} = 0.1 mm and 1 mm are also presented in Fig. 5 as dashed lines. The maximum pressure for particles with r < 1 mm is less than 3 Pa which is smaller than the corresponding tensile strength of the experimental particles. Thus, it is unlikely that a particle with r < 1 mm will disintegrate into smaller pieces because of pressure. For larger particles, the pressure can be higher. This happens as the sublimation creates a dry hot layer on the particle’s surface which helps to build up pressure at the sublimation front. As shown in Fig. 5, for the particle with r = 10 mm, the maximum pressure can reach the tensile strength of a particle made of aggregates with r_{a} = 0.1 mm. This indicates that the particleswith r > 10 mm may start fragmenting into smaller pieces but we expect the disintegration to stop at mmsized fragments. However, a detailed investigation of the fragmentation of particles is beyond the scope of this study and would require a detailed structural mechanical analysis. Albeit, we note that if such fragmentation would occur, it would have a significant effect on the sublimation rate as the insulating layer would be removed from time to time, allowing heat to directly access to the icy part.
Since the particles hold temperature gradients, outgassing of water vapour is anisotropic. This introduces a rocket force that accelerates the particles. The acceleration due to the rocket force is plotted in Fig. 6. The acceleration decreases with increasing particle size. This is clear as the sublimation rate is proportional to the illuminated area and the mass is proportionalto the volume, the acceleration roughly scales linearly with the inverse of size. For small particles, the acceleration is high but it only occurs for a few seconds, whereas for large particles, the acceleration is small but it acts for a considerabletime, which may have an effect on the particle’s trajectory. The direction of the acceleration depends on the rotationstate of the particle. For a slowly rotating particle, the force points towards the antisolar direction whereas for a fastrotating particle the direction is shifted towards the direction perpendicular to the rotation axis, e_{rot}, and to the axis towards the Sun, e_{sun}, as presentedin Fig. 7.
Finally, we studied how the surfacetovolume ratio, S, in Eq. (5) affects the results. Since sublimation requires that P_{sat} − P > 0 and it depends on the surfacetovolume ratio of the ice surface, S_{ice}, whereas condensation requires that P_{sat} − P < 0 and it depends on the total surfacetovolume ratio S_{tot} = S_{ice} + S_{dust}, we define (27)
If ice fully covers the dust grains, it is clear that S_{ice} = S_{tot}. If the particle is an aggregate of ice and dust grains, then S_{ice} < S_{tot}. For equisized icegrains: (28)
where v_{ice} is the volumetric filling factor of ice grains of radius r_{ice}. Thus, S depends on how ice is distributed inside the particle. We assumed that v_{ice} = 0.05 and compared four different ice distributions. We also assumed that the ice distribution does not change the scattering, absorption and emission properties. In the first case, ice covered all the dust grains giving S_{ice} = 1 × 10^{7} m^{−1}. In the second, the ice grains of radius r_{ice} = 0.2 μm were evenly distributed inside the particle giving S_{ice} = 7.5 × 10^{5} m^{−1}. In this case, the ice grains were approximately the same size as the dust grains. In the third, the ice grains were bigger, r_{ice} = 2 μm, than the dust grains and S_{ice} = 7.5 × 10^{4} m^{−1}. In the last case, the ice grains were much bigger than dust grains and r_{ice} = 20 μm.
Comparisons for the total ice volume fraction and the maximum pressure as a function of time are presented in Fig. 8. We observed that when S_{ice} is large enough, the system is saturated and increasing S_{ice} does not affect the results. This is because gas diffusion gives the upper limit to the total sublimation rate. When S_{ice} is small enough, the total sublimation rate is smaller and the system is no longer saturated and limited by gas diffusion but is limited by the ice surface area. In such a case, the received energy is used more to increase the temperature of the particle which in turn allows higher local pressure. In fact, as seen from Eq. (22), the thickness of the sublimation front depends on the length scale, h, and the diffusion coefficient, D, and is proportional to the factor, hSD^{−1}. When the factor is large, the sublimation front is thin compared to the particle’s size and the total sublimation rate depends on the macroscopic icy surface area rather than the microscopic surfacetovolume ratio, S. It is, thus, difficult to retrieve information on how the ice grains are mixed with the dust grains in the microscopic level from the total gas production rate.
Fig. 4 Ice volume fraction as a function of time for different particle sizes. 
Fig. 5 Maximum pressure as a function of time for different particle sizes. 
Fig. 6 Acceleration due to outgassing as a function of time for different particle sizes. 
Fig. 7 Acceleration vector components for a particle with r = 1 cm rotating 1 rpm about the axis, e_{rot}, perpendicular to the solar direction, e_{sun}, due to outgassing as a function of time. 
Fig. 8 Ice volume fraction and the maximum pressure as a function of time computed for varying ice surfacetovolume ratios, S_{ice}, corresponding to fully coated dust grains, and ice grains of sizes 0.2 μm, 2 μm, and 20 μm uniformly mixed with dust grains. 
6 Conclusions
We presented a novel thermophysical model and its numerical solution for icy cometary dust particles. The model employs the radiative transfer with reciprocal transactions for the radiative heat transport and the Boltzmann equation together with the HertzKnudsen formula for gas transport solved by the DSMC methods. Energy changes related to absorption, thermal emission, and phase changes were incorporated as source terms into the energy balance equation with the Fourier’s heat conduction and solved by the FEM.
The developed method allows for thermal analysis of up to cmsized particles composed of submicrometresized grains mixed with ice with a moderate computing power. The method can find applications in understanding the thermal physics of cometary dust particles and explaining experimental laboratory measurements, and help us to develop more accurate approximate thermal modelling methods.
We also showed that water ice sublimation inside large dust particles may generate enough pressure to reach the tensile strength and, thus, to possibly disintegrate the particles when comets are close to the Sun. Also, outgassing of water vapour can play a crucial role in the dynamics of icy dust particles after having been lifted off the comet’s nucleus. Both of these processes should leave observable effects on the remote observables via changes in the size and spatial distributions.
Finally, we showed that the thermal evolution of large dust particles is quite insensitive to the microscopic mixing scale of dust and ice, which effectively makes modelling easier. On the other hand, getting information on the mixing scale and ratio of dust and icefrom the total gas production rate is not trivial.
Acknowledgements
This work has been funded by the ERC Starting Grant No. 757390 Comet and Asteroid ReShaping through Activity (CAstRA). Computational resources have been provided by Gesellschaft für Wissenschaftliche Datenverarbeitung mbH Göttingen (GWDG).
References
 Agarwal, J., A’Hearn, M. F., Vincent, J.B., et al. 2016, MNRAS, 462, S78 [CrossRef] [Google Scholar]
 Ahmadian, M. H., Roohi, E., Teymourtash, A., & Stefanov, S. 2019, Phys. Fluids, 31, 062007 [CrossRef] [Google Scholar]
 Blum, J., Gundlach, B., Mühle, S., & TrigoRodriguez, J. 2014, Icarus, 235, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Brisset, J., Heißelmann, D., Kothe, S., Weidling, R., & Blum, J. 2016, A&A, 593, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503 [Google Scholar]
 Fanale, F. P., & Salvail, J. R. 1984, Icarus, 60, 476 [NASA ADS] [CrossRef] [Google Scholar]
 Gicquel, A., BockeléeMorvan, D., Zakharov, V. V., et al. 2012, A&A, 542, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gicquel, A., Vincent, J.B., Agarwal, J., et al. 2016, MNRAS, 462, S57 [CrossRef] [Google Scholar]
 Gundlach, B., Fulle, M., & Blum, J. 2020, MNRAS, 493, 3690 [CrossRef] [Google Scholar]
 Huebner, W., & Markiewicz, W. 2000, Icarus, 148, 594 [NASA ADS] [CrossRef] [Google Scholar]
 Huebner, W. F., Benkhoff, J., Capria, M.T., et al., eds. 2006, Heat and Gas Diffusion in Comet Nuclei [Google Scholar]
 Jäger, C., Mutschke, H., & Henning, T. 1998, A&A, 332, 291 [NASA ADS] [Google Scholar]
 Kelley, M. S., Lindler, D. J., Bodewits, D., et al. 2013, Icarus, 222, 634 [NASA ADS] [CrossRef] [Google Scholar]
 Kelley, M. S., Lindler, D. J., Bodewits, D., et al. 2015, Icarus, 262, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Lichtenegger, H., & Kömle, N. 1991, Icarus, 90, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Markkanen, J., & Agarwal, J. 2019, A&A, 631, A164 [CrossRef] [EDP Sciences] [Google Scholar]
 Markkanen, J., Penttilä, A., Peltoniemi, J., & Muinonen, K. 2015, Planet. Space Sci., 118, 164 [CrossRef] [Google Scholar]
 Markkanen, J., Agarwal, J., Väisänen, T., Penttilä, A., & Muinonen, K. 2018, ApJ, 868, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Prialnik, D. 2004, Comets II, 1, 359 [Google Scholar]
 Skorov, Y., & Blum, J. 2012, Icarus, 221, 1 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Crosscut of an example computational mesh generated by the twolevel Voronoi partitioning. The yellow tetrahedra correspond to the particle with a specified microstructure and the blue ones correspond to free space. 

In the text 
Fig. 2 Pressure as a function of distance from the centre of the sphere calculated by the DSMC and FEM. The DSMC 2 solution has twice the number of the DSMC molecules than the DSMC 1 solution to demonstrate convergence. 

In the text 
Fig. 3 Comparisons between thermal modelling results obtained by using the DSMC and FEM solutions for gas transport. The test particle has an irregular shape, with the volume equivalent radius of 0.5 mm. The surfacetovolume ratio was set to S = 1 m^{−1}. 

In the text 
Fig. 4 Ice volume fraction as a function of time for different particle sizes. 

In the text 
Fig. 5 Maximum pressure as a function of time for different particle sizes. 

In the text 
Fig. 6 Acceleration due to outgassing as a function of time for different particle sizes. 

In the text 
Fig. 7 Acceleration vector components for a particle with r = 1 cm rotating 1 rpm about the axis, e_{rot}, perpendicular to the solar direction, e_{sun}, due to outgassing as a function of time. 

In the text 
Fig. 8 Ice volume fraction and the maximum pressure as a function of time computed for varying ice surfacetovolume ratios, S_{ice}, corresponding to fully coated dust grains, and ice grains of sizes 0.2 μm, 2 μm, and 20 μm uniformly mixed with dust grains. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.