Thermophysical model for icy cometary dust particles

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 ice-free dust particles to ice-containing 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 centimeters 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 mm-sized fragments due to internal pressure build-up. Particles with r<1 cm lose their ice content within minutes. Hence, we expect that only particles with r>1cm may demonstrate sustained sublimation and the resulting outgassing forces.


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, ice-free particles (Gundlach et al. 2020). In addition, other mechanisms, such as the crystallisation of amorphous ices and impacts, may be capable of ejecting smaller icy particles. The ice-containing particles ejected from the surface may behave like mini-comets 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(Kelley et al. , 2015 and 67P/Churyumov-Gerasimenko . 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 wavelength-independent surface absorption and emission of electromagnetic radiation, the radiative heat transfer is treated as a one-dimensional 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 signifi-cant errors. The thermal evolution of small cometary dust particles is often modelled by assuming isothermal spherical particles and using the analytical Lorenz-Mie 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. 2012Gicquel et al. , 2016. This is a valid assumption for small compact particles, but inside porous particles larger than tens of micrometers, 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 analyzing the thermal properties of ice-free 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 ice-containing 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 finite-element 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 organised as follows. In Section 2, we present the governing equations describing the physics of the problem. The numerical methods we used are explained in Section 3. Section 4 introduces the particle model. In Section 5, we present the validation of the developed method against the continuous gas flow model and apply the method to study thermal properties of A&A proofs: manuscript no. AA39092 icy cometary dust particles at 1.35 AU. Finally, we present our conclusions in Section 6.

Governing equations
Recently, we introduced a thermal model for the aggregated icefree particles composed of submicrometer-sized dust grains in Markkanen & Agarwal (2019). Here, we will extend the method for dust particles containing ice inclusions. We assume that the dust and ice grains 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.

Energy balance
The energy balance equation for the solid phase is written as where c p is the specific heat capacity, ρ s is the density, T is the temperature, and κ is the conductive heat transfer coefficient. The right-hand 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 where q is the sublimation rate and H is the latent heat. By defining the saturation vapour pressure as and using the experimentally obtained coefficients, A and B (Fanale & Salvail 1984), the heat of sublimation is obtained via the Clausius-Clapeyron relation and given by 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 where S is the surface-to-volume ratio in the microscale, and P is the pressure in the micropores.

Mass balance
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: 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 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 and the gas flux as We apply the DSMC method to solve Equation (8), as described in Section 3.

Numerical solution
We discretise the domain of interest Ω with tetrahedral elements. Then we employ the FEM for the spatial dimension of the energy balance equation (1) with the nodal testing w n and basis functions u m . For the temporal dimension we will use the central finite difference formula and interpolate the other terms at t + 1 2 as This results into the Crank-Nickolson scheme for the unknown coefficient vector x t+1 given by where the mass and the stiffness matrices are defined as and the force vector as Equation (14) is strongly non-linear because of the force term. Hence, we use an iterative method with under-relaxation 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 molecule-molecule collisions is much larger than the mean free path of molecule-grain collisions. The grains are assumed to be much heavier than the gas molecules and, thus, they are stationary, giving us the so-called dusty gas model. We consider the molecule-grain 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 solid-phase temperature-dependent Maxwell-Boltzmann distribution. Then the trajectories of the DSMC molecules are traced within a time step. To account for the collision operator, C, the molecule-grain 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 Maxwell-Boltzmann 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 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 non-linearity 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 where m t is the ice mass at time t.

Particle model
We used a particle microstructure model derived by fitting the phase function of the coma of 67P/Churyumov-Gerasimenko 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 submicrometer-sized organic monomers and micrometer-sized 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 kgm −3 , are estimated from the porosity Φ and the silicate-to-organic 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, P rm 1 = 1, and P rm 2 = 0.4, 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 two level 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.

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 where D is the diffusion coefficient. Here, the gas density ρ g is defined as the gas mass per total volume including dust grains.
Using the continuity equation, we have Fig. 1. Crosscut of an example computational mesh generated by the two-level Voronoi partitioning. The yellow tetrahedra correspond to the particle with a specified microstructure and the blue ones correspond to free space. Writing the gas density in terms of pressure inside the micropores and using the ideal gas assumption, leads to the following equation for the gas density The diffusion coefficient for the random walk in three dimensions with the Maxwell-Boltzmann mean velocity v th is given by Furthermore, by assuming the half Maxwell-Boltzmann 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 Equation (22), supplemented with the boundary condition (24) can be solved with the FEM analogously to the FEM solution of the energy balance equation (14) by changing the coefficients and the right-hand 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 barycenter of each tetrahedron as a function of the distance from the center 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, Equation (22) is reduced to the Helmholtz equation in a steady state with the wavenumber squared given by The fundamental solution, that is, the Green's function for the Helmholtz equation is written as 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 S D −1 is large. This also means that in order to solve Equation (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 irregularly-shaped 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 surfaceto-volume 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 1rpm 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.

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 contain five volume percent 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 mm-sized particles. We observed that the particles with 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 non-linearity of Equation (14). A very small time step τ is required to solve the non-linear equation prohibiting long time evolution simulations. Thus, using a more sophisticated non-linear 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 T eff ≈ (r a /1 mm) −2/3 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 particles with r > 10 mm may start fragmenting into smaller pieces but we expect the disintegration to stop at mm-sized 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 proportional to 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 considerable time, which may have an effect on the particle's trajectory. The direction of the acceleration depends on the rotation state of the particle. For a slowly rotating particle, the force points towards the antisolar direction whereas for a fast rotating 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 presented in Fig. 7.
Finally, we studied how the surface-to-volume ratio, S , in Equation (5) affects the results. Since sublimation requires that P sat − P > 0 and it depends on the surface-to-volume ratio of the ice surface, S ice , whereas condensation requires that P sat − P < 0 and it depends on the total surface-to-volume ratio S tot = S ice + S dust , we define 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 ice grains: 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  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 surface-to-volume ratio was set to S=1 m −1 .  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 giv- ing 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 . Acceleration vector components for a particle with r = 1cm rotating 1 rpm about the axis, e rot , perpendicular to the solar direction, e sun , due to outgassing as a function of time.
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 Equation (22), the thickness of the sublimation front depends on the length scale, h, and the diffusion coefficient, D, and is proportional to the factor, hS D −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 surface-to-volume 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. S ice (20 m) = 7.5*10 3 m -1 Fig. 8. Ice volume fraction and the maximum pressure as a function of time computed for varying ice surface-to-volume 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.

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 Hertz-Knudsen 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 cm-sized particles composed of submicrometer-sized 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 A&A proofs: manuscript no. AA39092 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 ice from the total gas production rate is not trivial.