Thermal conductivity of porous aggregates
^{1} Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro, Tokyo 1528551, Japan
email: arakawa.s.ac@m.titech.ac.jp
^{2} Astronomical Institute, Tohoku University, Aoba, Sendai 9808578, Japan
^{3} Division of Theoretical Astronomy, National Astronomical Observatory of Japan, Mitaka, Tokyo 1818588, Japan
Received: 27 October 2017
Accepted: 16 November 2017
Context. The thermal conductivity of highly porous dust aggregates is a key parameter for many subjects in planetary science; it is not yet fully understood, however.
Aims. We investigate the thermal conductivity of fluffy dust aggregates with filling factors lower than 10^{1}.
Methods. We determined the temperature structure and heat flux of the porous dust aggregates calculated through Nbody simulations of static compression in the periodic boundary condition.
Results. We derive an empirical formula for the thermal conductivity through the solid network k_{sol} as a function of the filling factor of dust aggregates φ. The results reveal that k_{sol} is approximately proportional to φ^{2}, and the thermal conductivity through the solid network is significantly lower than previously assumed. In light of these findings, we must reconsider the thermal histories of small planetary bodies.
Key words: conduction / radiative transfer / comets: general / meteorites, meteors, meteoroids
© ESO, 2017
1. Introduction
Understanding the physical parameters of dust aggregates is important in planetary science. Specifically, the thermal conductivity of dust aggregates is key for determining the thermal evolution of planetary bodies, which influences the thermal evolution pathways of both rocky and icy planetesimals (e.g., Henke et al. 2012; Sirono 2017). The thermal evolution and activity of cometary nuclei also depend on the thermal conductivity of icy aggregates (e.g., Haruyama et al. 1993; GuilbertLepoutre & Jewitt 2011).
The thermal conductivity of dust aggregates depends on many parameters, and many previous experimental studies have researched the thermal conductivity of dust aggregates with filling factors above 10^{1}. The thermal conductivity of porous aggregates in vacuum is given by two terms: the thermal conductivity through the solid network k_{sol}, and the thermal conductivity owing to radiative transfer k_{rad}. Krause et al. (2011) showed that the thermal conductivity through the solid network k_{sol} is exponentially dependent on the filling factor of dust aggregates φ for 0.15 <φ < 0.54, and concluded that the coordination number of monomer grains C influences the efficiency of heat flux within the aggregates. Sakatani et al. (2016) revealed that k_{sol} is also dependent on the contact radius between monomers r_{c}. The thermal conductivity owing to radiative transfer k_{rad} is affected by the temperature of dust aggregates T and the mean free path of photons l_{p} (e.g., Schotte 1960; Merrill 1969). Moreover, l_{p} depends on R and φ when we apply the geometrical optics approximation for the evaluation of l_{p} (e.g., Skorov et al. 2011; Gundlach & Blum 2012).
There are also several theoretical studies on the thermal conductivity of dust aggregates (e.g., Chan & Tien 1973; Sirono 2014; Sakatani et al. 2017). However, no previous research has been conducted on the thermal conductivity of porous aggregates with filling factors lower than 10^{1}, although Kataoka et al. (2013b) and Arakawa & Nakamoto (2016) revealed that the collisional growth of dust aggregates leads to planetesimal formation via highly porous aggregates with filling factors much lower than 10^{1}. Therefore, the purpose of this study is to investigate the thermal conductivity and thermal evolution of fluffy dust aggregates in protoplanetary disks.
In this Letter, we calculate thermal conductivity through the solid network k_{sol} for highly porous aggregates with filling factors in the range of 10^{2} to 10^{1}. We use the snapshot data of Kataoka et al. (2013a) to calculate k_{sol}. We then validate our results through a comparison with the experimental data of Krause et al. (2011). We also derive the thermal conductivity owing to radiative transfer k_{rad} for porous aggregates of submicronsized monomers. Our results show that the thermal conductivity of highly porous aggregates is significantly lower than previously assumed.
2. Method
2.1. Arrangement of monomer grains
The arrangement of monomer grains depends on the coagulation history of the aggregates. During initial dust aggregate coagulation in protoplanetary disks, both experimental (e.g., Wurm & Blum 1998) and theoretical (e.g., Kempf et al. 1999) studies have shown that hitandstick collisions lead to the formation of fractal aggregates with a fractal dimension D ~ 2, which is called ballistic clustercluster aggregation (BCCA; Meakin 1991). Furthermore, Kataoka et al. (2013a) performed threedimensional numerical simulations of static compression of BCCA aggregates consisting of 16 384 spherical grains using a periodic boundary condition. In this study, we use snapshots of the compressed BCCA aggregates calculated by Kataoka et al. (2013a).
2.2. Temperature structure of the dust aggregate
To calculate the thermal conductivity through the solid network of an aggregate k_{sol}, we have to determine the temperature of each grain in a cubic periodic boundary. We calculated the temperature of each grain using the method of Sirono (2014). Here, we considered onedimensional heat flow from the lower boundary plane to the upper boundary plane. There are three choices regarding the pair of lower and upper planes, and we calculated k_{sol} from three directions. Then, we averaged these values for each snapshot.
We defined R as the monomer radius and L^{3} as the volume of each cubic space. The location of the ith grain (x_{i},y_{i},z_{i}) satisfies  x_{i}  <L/ 2,  y_{i}  <L/ 2, and  z_{i}  <L/ 2 for i = 1,2,...,N, where N = 16 384 is the number of grains in the periodic boundary. A sketch of a dust aggregate in a cubic periodic boundary is shown in Fig. 1. Here, we assumed that heat flow occurs along the zdirection. The grains located in −L/ 2 <z_{i} < −(L/ 2−R) are on the lower boundary (number 1 in Fig. 1), and the grains located in + (L/ 2−R) <z_{i} < + L/ 2 are on the upper boundary (number 40 in Fig. 1). When the ith grain was located on the lower (upper) boundary, we added a new grain on the upper (lower) boundary. The location of the new grain is (x_{i},y_{i},z_{i} + L) if the ith grain is located on the lower boundary (location X in Fig. 1) and (x_{i},y_{i},z_{i}−L) if the ith grain is located on the upper boundary (location Y in Fig. 1). We set the temperature of grains located on the lower (number 1 and location Y in Fig. 1) and upper (number 40 and location X in Fig. 1) boundary as T_{0} + ΔT/ 2 and T_{0}−ΔT/ 2, respectively.
Fig. 1 Sketch of a dust aggregate in a cubic periodic boundary. The temperature of grains located on the lower (number 1 and location Y) and upper (number 40 and location X) boundary is set to T_{0} + ΔT/ 2 and T_{0}−ΔT/ 2, respectively. The temperature of each grain is calculated by solving Eq. (1) simultaneously for each grain. 

Open with DEXTER 
Heat flows through the monomermonomer contacts, and for steadystate conditions, the equation of heat balance at the ith grain is given by (1)where F_{i,j} is the heat flow from the jth grain to the ith grain, given by (2)where H_{c} is the heat conductance at the contact of two grains, and T_{i} and T_{j} are the temperatures of the ith and jth grains, respectively. We considered the contacts not only inside the periodic boundary but also on the side boundaries (e.g., the contacts between numbers 9 and 10 and between numbers 21 and 22 in Fig. 1). The heat conductance at the contact of two grains H_{c} is (Cooper et al. 1969) (3)where k_{mat} is the material thermal conductivity and r_{c} is the contact radius of monomer grains. The contact radius r_{c} depends on the monomer radius R and the material parameters (Johnson et al. 1971). The heat conductance within a grain H_{g} is also given by (Sakatani et al. 2017) (4)However, we neglected the effect of H_{g} because H_{g} is sufficiently higher than H_{c} for (sub)micronsized grains. Therefore, the temperature structure of the aggregate in the cubic periodic boundary can be calculated by solving Eq. (1) simultaneously for all N grains, except for lower and upper boundary grains.
2.3. Thermal conductivity through the solid network
After we obtained the temperature structure, we calculated the total heat flow at the upper boundary ∑ _{upper}F_{i,j}, where we took the sum of contacts between the upper boundary ith grain and internal jth grain (for the case of Fig. 1, ∑ _{upper}F_{i,j} = F_{X,27} + F_{40,39}). The total heat flow at the upper boundary ∑ _{upper}F_{i,j} can be rewritten using the thermal conductivity through the solid network k_{sol} as (5)In this study, we discuss k_{sol} as a function of the filling factor φ, and rewrite L using φ as (6)Therefore, we obtain k_{sol} as a function of φ as follows: (7)where f(φ) is a dimensionless function of φ. Note that the total heat flow at the lower boundary −∑ _{lower}F_{i,j} is clearly equal to ∑ _{upper}F_{i,j} considering the heat balance.
3. Results
Here, we present the dimensionless function f(φ) for nine snapshots from three runs and three densities (Table 1). We calculated f(φ) in three directions for each snapshot and took the arithmetic mean values. We note that the compressed BCCA aggregates might be isotropic if the number of monomer grains is sufficiently large; thus, we only discuss the mean values of f(φ).
Fig. 2 Fitting of the dimensionless function of thermal conductivity f(φ) as a function of the filling factor φ. The green dashed line is the bestfit line, and the magenta solid line represents the simple function f(φ) = φ^{2}. 

Open with DEXTER 
Results of the numerical calculation.
Figure 2 shows the dimensionless function f(φ) as a function of the filling factor φ. The bestfit line given by the leastsquares method (green dashed line) is (8)Hereafter, we use the following more simple relationship between f(φ) and φ (magenta solid line) (9)Sakatani et al. (2017) predicted that f(φ) = (2 /π^{2})Cφ, where C is the coordination number. For highly porous aggregates, the coordination number C is approximately two, and the filling factor dependence on C is weak. Hence, f(φ) would be proportional to φ in the model of Sakatani et al. (2017); however, in reality, f(φ) is approximately proportional to φ^{2}.
In the context of the thermal conductivity of colloidal nanofluid and nanocomposites, Evans et al. (2008) revealed that thermal conductivity is strongly affected by the fraction of linear chains that contribute to heat flow in the aggregates. The contributing grains are called backbone grains, and noncontributing grains are called deadend grains (numbers 8, 32, and 33 in Fig. 1; Shih et al. 1990). We will discuss the effects of different fractions of backbone and deadend grains on the thermal conductivity in future research.
By comparing our model to the experimental data of Krause et al. (2011), we can confirm the validity of our model (Fig. 3). The magenta line in Fig. 3 represents the calculated thermal conductivity from Eqs. (7) and (9), the blue curve is the exponential fitting of experimental data, k_{sol} = 1.4e^{7.91(φ−1)} W m^{1} K^{1} (Krause et al. 2011), and the black dashed line is a model commonly used to study the thermal evolution of planetary bodies, that is, k_{sol} = φk_{mat} (e.g., Sirono 2017). Both experimental (crosses) and numerical data (squares, triangles, and circles) are plotted.
Fig. 3 Experimental (crosses) and numerical data (squares, triangles, and circles) of thermal conductivity through the solid network k_{sol}. Our model (magenta line) was compared to the experimental fitting data of Krause et al. (2011, blue curve). 

Open with DEXTER 
When we consider the dust aggregates of (sub)micronsized monomers, the contact radius between monomers r_{c} is given by (Johnson et al. 1971) (10)where γ = 25 mJ m^{2}, ν = 0.17, and Y = 54 GPa are the surface energy, Poisson’s ratio, and Young’s modulus of SiO_{2} grains (Wada et al. 2007). We set R = 0.75 μm and k_{mat} = 1.4 W m^{1} K^{1} to the same values as Krause et al. (2011). Figure 3 clearly shows that our empirical model is applicable to the k_{sol} of porous aggregates with filling factors of φ ~ 0.1. Moreover, our model is applicable not only for φ ≲ 0.1, but also for the range 0.1 ≲ φ ≲ 0.5.
We note that most of the experimental data of k_{sol} were fitted using exponential functions of φ, which pass the material thermal conductivities (e.g., Krause et al. 2011; Henke et al. 2013). However, when we consider the thermal conductivity through the solid network, that is, thermal conductivity limited by the necks between two monomers, k_{sol} must be given by k_{sol} ~ (r_{c}/R)k_{mat}, even for dense dust aggregates whose filling factors are close to unity (e.g., Chan & Tien 1973).
4. Discussion
Finally, we evaluated the total thermal conductivity of porous icy aggregates under vacuum conditions. The thermal conductivity through the solid network k_{sol} is given by (11)and thermal conductivity owing to radiative transfer k_{rad} is given by (Merrill 1969) (12)where σ_{SB} = 5.67 × 10^{8} W m^{2} K^{4} is the StefanBoltzmann constant. We calculated the mean free path of photons in fluffy aggregates of submicronsized grains l_{p} as follows: (13)where κ_{abs} and κ_{sca} are the absorption and scattering mass opacities of monomers, respectively, and ρ_{mat} = 1.68 × 10^{3} kg m^{3} is the material density. Here, the composition of icy dust aggregates is consistent with Pollack et al. (1994). The total mass opacity of submicronsized monomers κ_{abs} + κ_{sca} is hardly dependent on the wavelength of the thermal radiation λ = 2.9 × 10^{3} (T/ K)^{1} m for 10^{6} m <λ < 10^{4} m, and κ_{abs} + κ_{sca} is on the order of 10^{2} m^{2} kg^{1} (e.g., Kataoka et al. 2014). Then we set l_{p} = 10^{5}φ^{1} m in this study. We note that for the case of fluffy aggregates of submicronsized monomers, the wavelength λ is larger than the monomer radius R, even if the temperature is on the order of 10^{3} K. Hence, we cannot apply the geometrical optics approximation to evaluate l_{p}.
Figure 4 shows the k_{sol} of crystalline and amorphous icy aggregates, k_{sol,cr} and k_{sol,am}, and the thermal conductivity owing to radiative transfer k_{rad} for aggregates composed of icy monomers with a radius of R = 0.1 μm and temperature of T = 40 K. We set γ = 100 mJ m^{2}, ν = 0.25, and Y = 7 GPa for icy grains (Wada et al. 2007). The material thermal conductivities of crystalline and amorphous grains, k_{mat,cr} and k_{mat,am}, are given by k_{mat,cr} = 5.67 × 10^{2} (T/ K)^{1} W m^{1} K^{1} and k_{mat,am} = 7.1 × 10^{8} (T/ K) W m^{1} K^{1}, respectively (Haruyama et al. 1993).
Fig. 4 Comparison between k_{sol} (magenta solid line for crystalline icy aggregates and magenta dashed line for amorphous) and k_{rad} (blue dashed line). The monomer radius is R = 0.1 μm, and the temperature is T = 40 K. 

Open with DEXTER 
For the case of φ < 4 × 10^{3}, the thermal conductivity owing to radiative transfer k_{rad} is higher than the thermal conductivity through the solid network of crystalline icy aggregates k_{sol,cr}, even when the temperature is sufficiently low (T = 40 K). If the total thermal conductivity k_{sol} + k_{rad} does not change substantially when crystallization occurs, then the internal temperature of icy planetesimals could still increase after crystallization, which might cause runaway crystallization due to latent heat. In addition, when the temperature of an icy planetesimal increases, sintering can proceed inside the icy aggregate before monomer grains evaporate or melt (e.g., Sirono 2017). We cannot evaluate the contact radius r_{c} from Eq. (10) when aggregates are sintered, and k_{sol} increases linearly as a consequence of the increase of r_{c}. Sintering might also affect the mechanical strength of aggregates (e.g., Omura & Nakamura 2017) and the critical velocity for collisional growth (e.g., Sirono & Ueno 2017). Therefore, the growth pathways of icy planetesimals might be altered by sintering of icy aggregates. Not only icy planetesimals, but also rocky aggregates could experience sintering before growing into dmsized bodies, which might explain the retainment of chondrules inside fluffy aggregates (Arakawa 2017). We also note that the total thermal conductivity might be controlled by the thermal conductivity due to gas diffusion for the case of fluffy aggregates in high gas density environments (e.g., the innermost region of protoplanetary disks and/or planetary surfaces; Piqueux & Christensen 2009).
In conclusion, we have revealed the filling factor dependence of the thermal conductivity of porous aggregates. We showed that the thermal conductivity of highly porous aggregates is significantly lower than previously assumed. In future work, we will reexamine the growth pathways of planetesimals in protoplanetary disks, and combine this with a density and thermal evolution analysis.
Acknowledgments
We thank the referee Gerhard Wurm for constructive comments. This work is supported by JSPS KAKENHI Grant (15K05266). S.A. is supported by the GrantinAid for JSPS Research Fellow (17J06861).
References
 Arakawa, S. 2017, ApJ, 846, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Arakawa, S., & Nakamoto, T. 2016, ApJ, 832, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, C. K., & Tien, C. L. 1973, Journal of Heat Transfer, 95, 302 [CrossRef] [Google Scholar]
 Cooper, M. G., Mikic, B. B., & Yovanovich, M. M. 1969, International Journal of Heat and Mass Transfer, 12, 279 [CrossRef] [Google Scholar]
 Evans, W., Prasher, R., Fish, J., et al. 2008, International Journal of Heat and Mass Transfer, 51, 1431 [CrossRef] [Google Scholar]
 GuilbertLepoutre, A., & Jewitt, D. 2011, ApJ, 743, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Gundlach, B., & Blum, J. 2012, Icarus, 219, 618 [NASA ADS] [CrossRef] [Google Scholar]
 Haruyama, J., Yamamoto, T., Mizutani, H., & Greenberg, J. M. 1993, J. Geophys. Res., 98, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Henke, S., Gail, H.P., Trieloff, M., Schwarz, W. H., & Kleine, T. 2012, A&A, 537, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henke, S., Gail, H.P., Trieloff, M., & Schwarz, W. H. 2013, Icarus, 226, 212 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, K. L., Kendall, K., & Roberts, A. D. 1971, Proc. Roy. Soc. London Ser. A, 324, 301 [NASA ADS] [CrossRef] [Google Scholar]
 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013a, A&A, 554, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013b, A&A, 557, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kataoka, A., Okuzumi, S., Tanaka, H., & Nomura, H. 2014, A&A, 568, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kempf, S., Pfalzner, S., & Henning, T. K. 1999, Icarus, 141, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, M., Blum, J., Skorov, Y. V., & Trieloff, M. 2011, Icarus, 214, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Meakin, P. 1991, Rev. Geophys., 29, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Merrill, R. B. 1969, Nasa Technical Note, D5063 [Google Scholar]
 Omura, T., & Nakamura, A. M. 2017, Planet. Space Sci., in press [arXiv:1708.08038] [Google Scholar]
 Piqueux, S., & Christensen, P. R. 2009, J. Geophys. Res., 114, E09005 [NASA ADS] [Google Scholar]
 Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Sakatani, N., Ogawa, K., Iijima, Y.i.,Arakawa, M., & Tanaka, S. 2016, Icarus, 267, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Sakatani, N., Ogawa, K., Iijima, Y.i., et al. 2017, AIP Adv., 7, 015310 [NASA ADS] [CrossRef] [Google Scholar]
 Schotte, W. 1960, AIChE Journal, 6, 63 [CrossRef] [Google Scholar]
 Shih, W.H., Shih, W. Y., Kim, S.I., Liu, J., & Aksay, I. A. 1990, Phys. Rev. A, 42, 4772 [NASA ADS] [CrossRef] [Google Scholar]
 Sirono, S.i. 2014, Meteor. Planet. Sci., 49, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Sirono, S.i. 2017, ApJ, 842, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Sirono, S.i., & Ueno, H. 2017, ApJ, 841, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Skorov, Y. V., van Lieshout, R., Blum, J., & Keller, H. U. 2011, Icarus, 212, 867 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Wurm, G., & Blum, J. 1998, Icarus, 132, 125 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Sketch of a dust aggregate in a cubic periodic boundary. The temperature of grains located on the lower (number 1 and location Y) and upper (number 40 and location X) boundary is set to T_{0} + ΔT/ 2 and T_{0}−ΔT/ 2, respectively. The temperature of each grain is calculated by solving Eq. (1) simultaneously for each grain. 

Open with DEXTER  
In the text 
Fig. 2 Fitting of the dimensionless function of thermal conductivity f(φ) as a function of the filling factor φ. The green dashed line is the bestfit line, and the magenta solid line represents the simple function f(φ) = φ^{2}. 

Open with DEXTER  
In the text 
Fig. 3 Experimental (crosses) and numerical data (squares, triangles, and circles) of thermal conductivity through the solid network k_{sol}. Our model (magenta line) was compared to the experimental fitting data of Krause et al. (2011, blue curve). 

Open with DEXTER  
In the text 
Fig. 4 Comparison between k_{sol} (magenta solid line for crystalline icy aggregates and magenta dashed line for amorphous) and k_{rad} (blue dashed line). The monomer radius is R = 0.1 μm, and the temperature is T = 40 K. 

Open with DEXTER  
In the text 