Scattering, absorption, and thermal emission by large cometary dust particles: Synoptic numerical solution

Context: Remote light scattering and thermal infrared observations provide clues about the physical properties of cometary and interplanetary dust particles. Identifying these properties will lead to a better understanding of the formation and evolution of the Solar System. Aims: We present a numerical solution for the radiative and conductive heat transport in a random particulate medium enclosed by an arbitrarily shaped surface. The method will be applied to study thermal properties of cometary dust particles. Methods: The recently introduced incoherent Monte Carlo radiative transfer method developed for scattering, absorption, and propagation of electromagnetic waves in dense discrete random media is extended for radiative heat transfer and thermal emission. The solution is coupled with the conductive Fourier transport equation that is solved with the finite-element method. Results: The proposed method allows the synoptic analysis of light scattering and thermal emission by large cometary dust particles consisting of submicrometer-sized grains. In particular, we show that these particles can sustain significant temperature gradients resulting in the superheating factor phase function observed for the coma of comet 67P/Churyumov-Gerasimenko.


Introduction
The physical properties of cosmic dust particles, which include size, shape, porosity, and composition, can be constrained from remote light scattering and thermal emission observations. Observations of scattered solar and thermally emitted light have mostly been interpreted separately using different particle models and numerical methods (Kolokolova et al. 2004;Kimura et al. 2016). This may lead to contradictions in the retrieved particle properties as the interpretation procedure is not self-consistent. Using the same numerical method and the same particle model for scattered solar and thermally emitted light (i.e., synoptic modeling) results in more reliable estimates for the physical properties of dust particles. Such modeling, however, introduces various numerical challenges due to the multiscale nature of the problem. The solution is needed from ultraviolet to thermal infrared. In addition, cometary dust particles contain multiscale structures. They are thought to be large aggregates consisting of submicrometer-sized grains (Güttler et al. 2019).
Numerically exact electromagnetic techniques allow for analyses of scattering and radiative heat transfer properties of small aggregated particles. Combining the standard electromagnetic methods with the statistical fluctuation driven electrodynamic solvers, the radiative heat transfer problem can be solved directly from Maxwell's equations (Jin et al. 2017). This approach is computationally heavy, and consequently it is not applicable to particles that are much larger than the optical wavelengths or that contain a large number of individual grains. For sufficiently large objects such as asteroids and comet nuclei, a typical approach is to localize absorption and emission sources on the surfaces of the object using a proper boundary condition. Then, by introducing the effective conductive and radiative heat transfer coefficient, heat transport inside the object can be modeled (Huebner et al. 2006;Delbo et al. 2015). This simple model is a valid approximation only for one-dimensional problems, and it introduces significant errors when the object has locally threedimensional features (on length scales on the order of the mean free path length). Thus, large objects with small details cannot be treated rigorously with the existing computational techniques, due to the large computational time and memory requirement or inadequate approximations.
It is also possible to use the radiative transfer equation (RTE) in combination with the Fourier heat transport equation to deal with heat transfer in participating media (Howell et al. 2010;Modest 2013). A large number of different solution techniques to solve the RTE exist, for example the discrete ordinates (Hardy et al. 2016), finite-volume (Coelho 2014), finite-element methods (Richling et al. 2001), and Monte Carlo ray tracing methods (Howell 1998). Nonetheless, the standard RTE is only applicable to sparse discrete random media where the grains are in each other's far zone. Recently, a numerical framework for electromagnetic scattering, absorption, and propagation in dense discrete random media was introduced which extends the applicability range of the standard RTE (Muinonen et al. 2018;Markkanen et al. 2018a;Väisänen et al. 2019). Furthermore, Markkanen et al. (2018b) augmented the method with the coherent field contribution to account for surface reflections and refractions allowing the scattering analysis of millimetersized cometary dust particles consisting of submicrometer-sized grains.
In this paper we extend the numerical method developed in Markkanen et al. (2018b) for the radiative heat transfer problems. The radiative heat transfer solution is coupled with the transient conductive heat transfer equation, which is solved by the finite-element method (FEM). The presented approach allows us to analyze the scattering and thermal properties of particles consisting of aggregates of small grains over a wide wavelength range. Finally, we show that large cometary dust particles can hold significant temperature gradients, which explains the superheating phase function of the coma of comet 67P/Churyumov-Gerasimenko (67P) measured by the Rosetta/VIRTIS (Visible and Infrared Thermal Imaging Spectrometer) instrument.

Methods
In this section, we present our numerical framework for solving radiative and conductive heat transfer in a particulate random medium enclosed by an arbitrarily shaped boundary. We neglect heat convection and assume that external radiation and thermal fluctuations are the only heat sources. We also assume that the radiative heat transport is much faster than the conductive heat transport. Thus, we consider the radiative part as instantaneous time harmonic radiation.

Energy equation
Let us consider a finite body Ω ∈ R 3 characterized by the position r dependent density ρ, specific heat capacity c p , and conductivity κ in free space. Our goal is to solve the time t dependent temperature T under the influence of the external radiation source. The energy balance equation in Ω for the combined conductive and radiative heat transfer reads as with the boundary condition n · κ∇T = n · q r on ∂Ω, where q r denotes the total radiative flux, ∂Ω is the boundary of Ω, and n is the inward pointing normal vector on ∂Ω.

Finite-element solution for the energy equation
To solve the energy balance Eqs. (1) and (2), we apply the standard finite-element method to discretize the spatial component of the energy equation. Taking the inner product with the testing functions w and integrating by parts, the weak formulation reads as follows: is valid for all w ∈ H 1 (Ω). Here H 1 (Ω) denotes the space of square integrable functions whose derivatives are also square integrable. We note that Eq. (3) automatically satisfies the Neumann boundary condition (2). Next we discretize the time derivative using the finitedifference formula where the subscript t indicates the time step and τ is its size. Substituting (4) into (3) and using the implicit backward Euler time integration scheme, we obtain We discretize the domain Ω with linear tetrahedral elements and expand the unknown temperature T into the linear nodal basis functions u m as T ≈ m x m u m , where x m are the unknown coefficients. Using Galerkin's method, in which the testing functions and basis functions are identical w m = u m , we can write the mass and stiffness matrices as and the force vector as Finally, we arrive at the expression for the temperature coefficient vector x given by The force vector is strongly nonlinear with respect to temperature (F ∼ T 4 ). Hence, estimating τF x t+1 with τF x t would require using a very small time step τ to reach a converging numerical solution. Therefore, we solve Eq. (9) for x t+1 by using an iterative method. In the next section we discuss how to compute the force vector F in a densly packed particulate random medium enclosed by an arbitrarily shaped boundary ∂Ω with the Monte Carlo radiative transfer approach.

Radiative transfer solution
To calculate the absorbed and thermally emitted power (i.e., the force vector F in Eq. (9)), we employ an approximate version of the recently introduced radiative transfer with the reciprocal transactions (R 2 T 2 ) method for electromagnetic scattering and propagation in dense discrete random media (Muinonen et al. 2018;Markkanen et al. 2018a;Väisänen et al. 2019). R 2 T 2 is based on the Monte Carlo ray tracing of the order-of-scattering interactions in which the fundamental scatterers are volume elements containing a large number of small grains. Incoherent electromagnetic interactions within the volume elements and the order-of-scattering interactions between the volume elements are computed via the exact electromagnetic solver, namely the fast superposition T-matrix method (FaSTMM; Markkanen & Yuffa 2017).
Here we use the far-field approximation to compute the order-of-scattering interactions. Yet, the interactions within the volume elements are computed exactly. We also neglect the reciprocal scattering diagrams (i.e., those responsible for the coherent backscattering effect). This allows us to include the coherent field contribution in terms of reflections and refractions on the surfaces of a particle using geometric optics (Markkanen et al. 2018b;Muinonen et al. 2019).
Computations proceed as follows. First, we compute scattering properties for volume elements, as described in Sect. 2.3.1. This step is done separately for each wavelength and for each A164, page 2 of 7 material region, but it does not depend on the shape and size of the entire particle. Second, the computed scattering properties from the first step are used as input parameters in the combined Monte Carlo geometric optics radiative transfer solver, as discussed in Sect. 2.3.2. This step is computationally independent of the first step for different particle shapes and sizes.

Scattering properties of volume elements
Input parameters for our radiative transfer algorithm are the incoherent scattering phase function M ic 11 , incoherent albedo ω, incoherent mean free path , coherent effective refractive index m eff , and absorption cross section of the ensemble averaged volume element C abs .
We compute these parameters for spherical ensembleaveraged volume elements that contain a large number of small spherical monomers. The radius of the volume element R 0 should be large enough to capture the statistics of the random medium. First, we generate a sample medium into a periodic box by randomly depositing monomers into the box until the target volumetric filling factor v is reached. Second, we draw a sample spherical volume element from the box, as explained in detail by Väisänen et al. (2019). Third, we apply FaSTMM to compute the Jones matrix S i of such a volume element i. Repeating the process N times (N = 512 is used in this work), we compute the coherent Jones scattering matrix as and the incoherent Jones scattering matrix for the ith state as Averaging the square of the incoherent Jones matrix over the ensemble of all N sample volume elements, we obtain the ensemble-averaged one-one element of the Mueller matrix, composed of the incoherent scattering phase function and the incoherent scattering cross section where S 2 is the surface of the unit sphere. The incoherent albedo is defined as where C abs is the ensemble-average absorption cross section, and the incoherent mean free path length as in which V is the volume of the volume element. Finally, we solve the coherent effective medium parameter m eff . This is done by matching the coherent scattering cross section C c sca to the C sca of the equal-sized sphere using the Mie solution. C c sca is defined as in (13), but the incoherent Mueller matrix is replaced with the coherent matrix.

Monte Carlo ray tracing
We trace the rays in the same tetrahedral mesh as is used in the FEM. This makes it trivial to couple the radiative transfer solution (i.e., ∇ · q r ) to the FEM solution. We decompose the divergence of the total flux into two parts ∇ · q r = q sol + q em in which q sol is the absorption rate due to the solar radiation, and q em is the absorption plus emission rate due to the thermally fluctuating sources.
To compute q sol , we launch N e ray rays for a given wavelength band ∆ λ from outside Ω. The power of each incident ray is given by where P(λ) is the flux density of the incident radiation and G is the geometric cross section of the particle. Each ray has a specific position r and propagation directionk depending on the source. For thermally emitted radiation, to compute q em , N i ray rays are launched from each tetrahedron C t with a random position and direction. The emitted power per ray for a given wavelength is given by where T c t is the average temperature of the tetrahedron C t and N d is the number density of the volume elements in the tetrahedron C t . The Planck function in a medium with the refractive index m written in terms of the free space wavelength λ 0 reads as where h is the Planck constant, c is the speed of light in vacuum, and k B is the Stefan-Boltzmann constant.
Once we have the initial position r, directionk, and power of the ray E, we start tracing. If the ray hits the boundary ∂Ω, it reflects and refracts according to Snel's law and the power is updated from Fresnel's coefficients which are calculated using the effective refractive index m eff . The position and direction of the reflected and refracted rays are updated.
If the ray is in Ω, the distance to the next scattering event is generated as where Φ is a uniform random number within [0, 1]. The scattering event happens if the ray does not cross the boundary ∂Ω.
Absorbed power E abs = E(1 − ω) in the scattering event is added into the total energy of the tetrahedron where the scattering event occurs, and the ray's power is updated (E = Eω). A new propagation directionk is drawn from the cumulative distribution function generated from the phase function of the ensemble averaged incoherent volume element M ic 11 . Finally, if the ray does not scatter and it does not cross the boundary ∂Ω, intensity is collected and added to the total scattering phase function. This process is repeated until the ray's power has decreased under the predefined threshold.
The radiative transfer computation gives us C t ∇ · q r dV for each tetrahedron, which is then used to compute the force vector F using (8). In addition, since the escaped rays are recorded, the total scattering and thermal emission phase functions and cross sections are obtained for each wavelength.

Numerical results
In this section, we present some numerical results. To validate the method, we compare our solution with the known approximate solution for small loosely packed grains computed by employing the effective medium approximation (EMA) and the Mie theory. Then, we apply the method to explain the superheating factor phase function measured for the coma of comet 67P.

Comparison to the effective medium approximation
Let us first consider a spherical medium consisting of small spherical silicate grains of radius r = 20 nm. The grains are randomly positioned in a spherical domain Ω with radius R. The volumetric filling factor is ν = 0.05. The wavelength dependent refractive index for magnesium-iron silicate, from Dorschner et al. (1995), is plotted in Fig. 1. Here we assume that the refractive index is independent of temperature in the studied temperature range.
The external radiation source is the Sun, which is assumed to radiate as a perfect blackbody of temperature T = 5777 K. The solar flux P in the given wavelength band (λ 1 − λ 2 ) at the distance ∆ = 1.3 AU from the Sun is given by where r is the radius of the Sun. Since the grain size is small, we can calculate a reference solution at the equilibrium temperature using the steady-state energy balance equation, EMA, and the Mie theory. To compute the absorption cross section C abs , we use the Maxwell-Garnett mixing rule (Sihvola 2000) and the Mie theory. For an isothermal spherical particle, the steady-state energy balance equation reads as To compare the proposed method to the steady-state analysis with EMA (21), we ran the transient simulation until the   Figure 2 shows the total absorbed solar radiation and thermally emitted power as a function of wavelength calculated by the proposed method and EMA for three different particle sizes R = 0.01, 0.1, and 1.0 mm. We observed excellent agreement between the two methods; however, we expect a decrease in the accuracy of our method as the particle size becomes smaller. This is evident as we used geometric optics to approximate the coherent reflections and refractions on the surface of a particle.
In the computations the volume element radius is r 0 = 0.5 µm and the time step is τ = 0.1 s, τ = 1.0 s, and τ = 10 s for 0.01 mm, 0.1 mm, and 1.0 mm particles, respectively. We also assumed high thermal conductivity κ = 100 W m −1 K −1 to remove possible temperature gradients inside the body. Although it is not a realistic assumption for real dust, it serves as a validation by simplifying the problem. It is worth noting that the density ρ = 1000 kg m −3 and the specific heat capacity c p = 500 J kg −1 K −1 do not affect the equilibrium temperature as ∂T ∂t = 0 in (1) for the steady-state solution; they only affect the time it takes to reach equilibrium. Figure 3 shows the temperature evolution of the particles when the initial temperature T 0 = 200 K. We see that the transient temperature reaches the equilibrium temperature computed by using the steady-state approach with EMA.
Next we studied larger grains with r = 100 nm and increased the filling factor to v = 0.3. These values are commonly used to model cometary dust (Kolokolova et al. 2004). It is evident from Figs. 4 and 5 that EMA leads to a different absorption rate, which consequently affects the equilibrium temperature and thermal emission; the EMA is no longer a valid approximation at shorter wavelengths as the grain size is approximately the same as the wavelength. Thus, care should be taken when applying EMAs to interpret infrared observations of cometary dust particles consisting of 100 nm-sized grains.

Application to the coma of 67P
We recently introduced a particle model that explains the scattering phase functions of the coma of comet 67P measured by the Rosetta/OSIRIS (Optical, Spectrocopic and Infrared Remote Imaging System) instrument (Markkanen et al. 2018b). The particle model consists of aggregated submicrometer-sized organic grains and micrometer-sized silicate grains. The shapes of the aggregates are Gaussian random spheres. It is thus interesting to study whether the same model can also explain the thermal infrared observations of the same comet. Since there is no refractive index available for cometary organic material over a wide wavelength range, we used the refractive index of carbonaceous dust analogues (Jäger et al. 1998). For silicates, we used the refractive index presented in Fig. 1 that corresponds to magnesium-iron silicate mineral (Dorschner et al. 1995).
We concentrate on the superheating factor phase function reported by Bockelée-Morvan et al. (2019) for the coma of the comet 67P. The superheating factor is defined as the ratio of the color temperature T c to the equilibrium temperature of an ideal blackbody T b = 278/∆ 2 [K] (i.e., S = T c /T b ). To obtain the color temperature T c for a given phase angle, Bockelée-Morvan et al. (2019) fitted the Planck function to the Rosetta/VIRTIS near-infrared spectra (3−5 µm). They found that the superheating factor phase function has a phase dependence indicating that the particles can hold temperature gradients if the particles are assumed to be randomly oriented.     Fig. 7. Difference between the maximum and minimum temperature inside the particle as a function of particle radius for an example particle. difference between the maximum and minimum temperature inside the particle with respect to the particle size is plotted in Fig. 7. It is evident that in order to model the measured superheating factor phase function slopes, the dust particles must sustain significant temperature gradients. Large temperature gradients, in turn, may have interesting effects on the dust dynamics in the inner coma via rocket forces caused by the asymmetrical sublimation of ices. Thermal break-down of particles may also play a role as the dust is transported from the inner to outer coma. Nevertheless, these interesting topics are beyond of the scope of this paper and are left for future research.
In our computations, the heat conduction coefficient is assumed to be κ = 0.00025 W m −1 K −1 , which corresponds to a realistic value for a porous dust aggregate (Krause et al. 2011;Arakawa et al. 2017;Sakatani et al. 2017). The heat conduction coefficient depends on porosity, grain size, and composition. Increasing the heat conduction coefficient flattens the superheating phase functions. Decreasing κ increases the slope until the radiative heat transport becomes more dominant. We also assumed that the angular velocities of the particles in the direction perpendicular to the Sun are slow (ω = 2π 50 0.1 mm R rad s −1 ),  Fig. 8. Modeled superheating factor phase function for different conductive heat transfer coefficients κ averaged over the differential powerlaw size distribution of index −3. The minimum cutoff is 10 µm and the maximum 100 µm.
hence the contribution of rotation to the temperature gradients is negligible. The slow rotation rate was introduced in order to speed up numerical convergence. Higher angular velocities decrease temperature gradients and flatten the superheating factor phase functions. However, it is important to note that if the rotation axis points towards the Sun, it does not affect the temperature gradients.
Finally, we averaged the modeled superheating phase functions over a differential power-law size distribution of index −3. The averaging results in a superheating factor phase function that is consistent with the measured value if the conductive heat coefficient κ = 0.00025 W m −1 K −1 , as demonstrated in Fig. 8. Such a low κ value suggests that the particles must be porous. The superheating factor phase function also indicates that the dust particles in the inner coma must be large. The dominating particle size range is around 10-100 µm. Smaller particles would show flatter superheating factor phase functions with higher absolute values whereas larger particles would result in a steeper slope and smaller absolute values. The size range is also consistent with other dust models explaining the intensity phase functions measured by Rosetta/OSIRIS (Moreno et al. 2018;Markkanen et al. 2018b).

Conclusions
We presented a synoptic numerical solution for scattering, absorption, and thermal emission by large aggregated dust particles consisting of submicrometer-sized grains. The numerical solution is self-consistent and allows for analysis of multiinstrument data sets. We solved the transient energy equation using the finite-element method in which the radiative heat transfer part is included as an additional forcing term. We developed a new algorithm, based on the Monte Carlo radiative transfer with reciprocal transactions framework, to compute the forcing term in the dense discrete random medium enclosed by an arbitrarily shaped surface. The method explicitly computes the radiative heat transfer part in contrast to the standard thermal models in cometary science where it is included in the heat conduction coefficient by assuming that heat transport is one dimensional. Thus, the developed method provides a full three-dimensional solution for radiative heat transfer problems. A164, page 6 of 7 We compared the solution of the proposed method to that obtained by employing the effective medium approximation and Mie theory. The solutions match if the grains are small and loosely packed, but deviate for grain sizes and packing densities that are typically assumed to form cometary dust particles. This implies that using EMAs to interpret thermal infrared observations of cometary dust may lead to unreliable results. Finally, we showed that particles 10−100 µm in size can reproduce the superheating factor phase function of the coma of comet 67P if the dust particles hold significant temperature gradients. This means that the heat conduction coefficient must be very low, hence particles must be porous. The presented analysis also suggests that the visible phase functions measured by Rosetta/OSIRIS and the superheating phase function from Rosetta/VIRTIS are consistent with each other.