Heat transfer in granular media with weakly interacting particles

We study the heat transfer in weakly interacting particle systems in vacuum. The particles have surface roughness with self-affine fractal properties, as expected for mineral particles produced by fracture, e.g., by crunching brittle materials in a mortar. We show that the propagating electromagnetic (EM) waves and the evanescent EM-waves, which occur outside of all solids, give the dominant heat transfer for large and small particles, respectively, while the phononic contribution from the area of real contact is negligible. As an application we discuss the heat transfer in rubble pile asteroids.

1 Introduction Granular materials can be described as homogeneous media in the continuum approximation, on length scales much larger then the particle sizes (diameter 2R).The effective thermal conductivity K is an important property of granular materials which may be very different from the thermal conductivity K 0 of a solid block made from the same material [1,2].Thus K depends on not only K 0 but on the size and shape of the particles, on the nature of the particle contact regions, and on the fraction of the total volume occupied by the particles, the so called filling factor (= 1 -macroporosity).It also depends on the environment such as temperature, gas pressure and humidity.
In this study we will assume that K << K 0 as is typically the case for granular solids in vacuum where there is no gas or fluid which could facilitate the heat transfer between the particles.In this case, if all the particles are small enough, the temperature in each particle may be nearly constant but the temperature change slightly from one particle to a nearby particle.The effective thermal conductivity is determined by the heat transfer between the particles and from dimensional arguments one expect K ≈ G R, where G is the thermal contact conductance relating the heat transfer rate between two particles to the temperature difference Q = G(T 0 − T 1 ).
Most studies of the heat transfer in granular materials have assumed spherical particles without surface roughness.However, all solids have surface roughness [3] which affect all contact mechanical [4], electrical [5] and thermal properties [6].In this paper we will study the influence of surface roughness on the heat transfer between particles in granular media.A similar study but using a very different formalism, and with less realistic types of surface roughness, was presented by Krüger et al (see Ref. [7]).
As an application we will consider the conduction of heat in asteroids [8,9] or in regolith, the loose material covering many surfaces of solar system bodies.Many asteroids consist of weakly interacting particles of different sizes (rubble pile asteroids).The effective thermal conductivity of a thin surface layer of asteroids can be obtained from the measured heat radiation, and are typically found to be a factor ∼ 0.01 times smaller than expected from a solid block made from the same material [10].We will show that the heat conduction is mainly due to the radiative (for large particles) and evanescent (for small particles) electromagnetic field, while the contribution from the area of real contact is negligible for weakly interacting particles.Our result constitutes a totally new, revised concept of what is conventionally called the "solid" thermal conductivity in granular matter studies.

Heat transfer in asteroids
All asteroids rotate and their surfaces are exposed to energy flux from the sun.Thus the asteroid surface temperature will vary in space and time and this can be studied using infrared detectors.As an example, Fig. 1(a) shows an optical image of the asteroid Benny and Fig. 1(b) the temperature distribution obtained from thermal emission.
On the average an asteroid must emit as much heat radiation as it absorb from the sun.Using this it is trivial to estimate the (average) asteroid surface temperature T which is (in equilibrium) also the nearly constant interior temperature.The total sun radiation power is Q ≈ 4 × 10 26 W. Most asteroids (radius r a ) are at distances r 0 ≈ (3 − 5) × 10 11 m from the sun.If we assume that all the radiation (photons) from the sun hitting an asteroid is absorbed by the asteroid then the absorbed power q = Qπr 2 a (4πr 2 0 ).The thermal radiation from the asteroids is the so-called "fast-rotator model" is σT 4 4πr 2 a , where σ ≈ 5.67×10 −8 W (m 2 K 4 ) is the Stefan-Boltzmann constant.Thus σT 4 4πr 2 a = Qπr 2 a (4πr 2 0 ) or , which gives T ≈ 200 K.The asteroids in the asteroid belt between Mars and Jupiter was formed at the same time as our solar system and is hence ≈ 4.6 × 10 9 years old.In Ref. [11] we showed that for a solid particle located in vacuum for ≈ 4.6×10 9 years at the temperature T ≈ 200 K all molecules bound to the particles by less than ≈ 1 eV will have desorbed.Hence we can assume that no mobile molecules exist which could form capillary bridges between the particles in asteroids.
The thermal emission of asteroids contains many important clues about their physical properties and the study of asteroid thermal emission (often referred to as thermal radiometry), together with optical brightness observations, is the primary source of known diameters and albedos (the measure of the diffuse reflection of solar radiation out of the total solar radiation), and the only established remote-sensing means of determining the crucial thermal inertia defined below.The principle of thermal radiometry is simple: Asteroids are heated up by absorption of sunlight, the absorbed energy is radiated off as thermal emission.This leads to a characteristic T (t, x) curve for each point x on the surface, or for the disc-averaged infrared flux for non-resolved observations.When a solid is exposed to a photon energy flux J(t) the temperature in the body will change with time.If the lateral variation of the energy flux J is slow, and the time variation fast, the lateral (parallel to the surface of the solid) heat diffusion can be neglected, and only heat diffusion normal to the surface (coordinate x) need to be considered.In this case the surface temperature depends only on the thermal inertia defined as Γ = (ρCK) 1 2 , where ρ is the bulk mass density, C the heat capacity and K the thermal conductivity.This is easily seen for the simplest case where the incident energy flux J(t) oscillate harmonically in time.Writing The boundary condition for x = 0: Thus the surface temperature depends only on the incident heat flux (amplitude J 0 and frequency ω) and on the thermal inertia Γ.During one oscillation in the energy flux (period 2π ω) the temperature penetrate a distance (the skin depth) l ≈ 1 Reγ ≈ K (Γ √ ω) into the solid.Asteroids or regolith consist of solid particles of different sizes (e.g., [12]) and the results presented above are only valid on length scales where the solid can be considered as effectively homogeneous, which would be length scales larger then the particle diameter for particles with the same size.Nevertheless, temperature map of the asteroids provide insight into their surface properties.Thus low thermal inertia is typically associated with layers of conducting wire resistor FIG. 3. Heat flow in a collection of spherical particles (radius R) containing a particle with bigger radius.The system can be replaced by an electric analogy consisting of perfect conduction wires (line segments) with resistors (black dots).
The big particles will reduce the flow resistance of the system by effectively reducing the number of resistors along some wires while keeping the number of wires (or flow channels) unchanged.Thus embedding big particles in a matrix of smaller particles may increasse the thermal conductivity, as indeed observed experimentally (see Fig. 8 in Ref. [1]).
dust, while high thermal inertia may indicate rocks on the surface.

Heat diffusion
We assume that the heat transfer between the particles is so slow that the temperature in each particle is approximately constant in space.We assume for simplicity that all the particles are spherical with equal radius R and forming a simple cubic lattice, and that the temperature depends only on the x-coordinate.Consider a particle at at position x at the temperature T (see Fig. 2).The heat transfer rate from the particle at x − ∆x into the particle at x is G[T (x − ∆x) − T (x)] and the heat flow rate out of the particle at x to the particle at where G is the (thermal) contact conductance.The net heat flow into the particle is This will give rise to a change in the temperature in the particle determined by where ρ is the particle mass density, C the heat capacity per unit mass and V = 4πR 3 3 the volume of the particle.Thus we get Using that ∆x = 2R we get the heat diffusion equation where the effective heat conductivity The prefactor γ = 3 π was derived for a simple cubic arrangements of the particles.Other arrangements of the spheres result in similar expressions for K, but with different prefactor γ (of order unity) (see Appendix A).We can also write where the heat diffusivity Note that the thermal diffusivity D and the effective heat conductivity K both have a factor 1 R. Thus, since for small particles with very rough surfaces G is nearly independent of the radius of the particles (typically, G ∼ R 0.2 , see Sec. 3), both D and K will increase as the size of the particles decrease.This counterintuitive result is easy to understand from the electric analogy shown in Fig. 3: decreasing the radius of the particles increases the number of contact points along the heat flux lines proportional to 1 R but at the same time the number of flux lines increases as 1 R 2 .As a result the heat resistivity will decrease as R 2 R = R and the heat conductivity increase as 1 R.For big particles the (radiative) black-body radiation will give a contribution to G proportional to R 2 and in this case the effective heat conductivity K will increase with the particle size as K ∼ R. For a system consisting of big particles surrounded by a matrix of small particles the effective heat conductivity may be determined by the small particles and the volume fraction they occupy in the mixture (see Fig. 3 for an electric analogy, and also Sec. 4).We The heat transfer between two particles in vacuum can occur via lattice vibrations (phonons) in the area of real contact (white arrows) or by electromagnetic radiation (photons) in the non-contact area (pink arrows).
note that the effective conductivity of a system with randomly distributed particles of different sizes and thermal conductivity could be studied using an effective medium approach.
4 The thermal contact conductance G Natural mineral particles, e.g., stone powder produced by crunching (involving fracture), are not perfect spherical but have very rough surfaces (see Fig. 4).This will drastically influence the heat contact conductance.In the most general case the heat transfer between two particles can occur via several different processes: (a) Contribution from the area of real contact.For insulators this corresponds to heat transfer via phonons (lattice vibrations) which can propagate from one solid to the other via the area of real contact.
(b) Heat radiation.Here one must in general consider both the propagating electromagnetic (EM) waves, which corresponds to the normal black-body radiation, and the evanescent EM-waves, which decay exponentially with the distance from the surfaces of the solids.The latter will dominate the heat transfer at short surface separation and will be very important for small particles.
(c) Heat transfer can occur in the surrounding atmosphere by heat conduction (or convection) in the gas.
(d) In the normal atmosphere fluid capillary bridges may form and heat diffusion in the fluid bridges will contribute to the heat transfer.
Here we are interested in the heat flow in asteroids and in these case there is no atmospheric gas and no capillary bridges so only processes (a) and (b) are relevant.

Area of real contact contribution to G
The thermal resistance of a contact is usually assumed to be due mainly to the constriction resistance.However, this assumes that the material at the interface interact as strongly as in the bulk.This may be the case in many practical (engineering) applications where the material in the contact regions are plastically deformed and where (for metals) cold welded regions form.However, here we are interested in the contact between particles which interact very weakly.In this case most of the contact resistance may be due to the weak coupling between the solids at the interface.We will now show that this is the case for the interaction between particles in asteroids.
For an atomic-sized contact the interfacial heat conductance at high enough temperature (see Appendix B and Ref. [6,[13][14][15]): where k B is the Boltzmann constant and η a damping due to phonon emission given by where c T is a transverse sound velocity, ρ the mass density of the solids and m the atomic mass.The parameter µ is given by (B6) (see also Ref. [13]) and depends only on the the Poisson ratio ν, or on the ratio c T c L between the transverse and longitudinal sound velocity.For c T c L = 1 2 one gets µ ≈ 0.13.The heat conductance G a for any temperature is obtained by replacing k B in (2) with C V 3 where C V is the heat capacity per atom.At high temperatures C V = 3k B (Dulong-Petit law) and we recover (2).Note that the heat capacity can be written as an energy fluctuation term (see Eq. (B7)) so the heat conductance (3) can be expressed as an energy fluctuation term times a damping (or inverse relaxation time) term η.
In (3) occurs the spring constant k = U ′′ (d 0 ) where U (d) is the interaction potential between surface groups on the two solids and d 0 the equilibrium separation.We assume a Lennard-Jones (LJ) interaction potential: Expanding in d − d 0 to second order gives Thus k = 72ǫ d 2 0 and using the same LJ-parameters as in Ref. [11], ǫ = 3.3 × 10 −21 J (or 0.02 eV) and d 0 = 0.3 nm, gives k ≈ 2.6 N m.Using the average atomic mass (the mass of SiO 2 divided by 3) m ≈ 3 × 10 −26 kg, the silica mass density ρ ≈ 2600 kg m 3 and the sound velocity c T ≈ 4000 m s we get η ≈ 10 11 s −1 .Thus the heat conductance G a = ηk B ≈ 10 −12 W K and the contact resistance R a = 1 G a ≈ 10 12 K W. We note that this thermal resistance is even larger then found for molecular junctions (where typically R ≈ 10 11 K W) involving small molecules between two gold electrodes, but in these cases the molecules are chemically attached to the metal electrodes [16,17].
The calculation above assumes implicitely that the bulk thermal conductivity is infinite.For finite thermal conductivity there will be an additional contact constriction resistance.We will show that the interfacial 5. The (phononic) interfacial contact resistance and the constriction resistance act in series and since the former is much larger the latter can be neglected.
resistance is much larger than the constriction resistance 1 (2r 0 K 0 ) where K 0 ≈ 1 W Km is the silica thermal conductivity and 2r 0 ≈ 1 nm the diameter of the contact area.Thus 1 (2r 0 K 0 ) ≈ 10 9 K W. Since the constriction resistance and the interfacial resistance R a act in series they add together.Hence the constriction resistance can be neglected and the contact conductance is determined by ( 2).This conclusion is summarized in Fig. 5.

Radiative contribution to G: theory
The heat flux per unit area between two black-bodies with flat surfaces (of area A 0 ∝ R 2 ) separated by d is given by the Stefan-Boltzmann law where T 1 and T 0 are the temperatures of solids 1 and 0, respectively, ̵ h the reduced Planck's constant and c the light velocity.Note that (4) is only valid if the surface separation d is larger than the wavelength λ of the emitted radiation.Since ck = ω where the wavenumber Hence the radiative heat transfer coefficient and the interfacial heat conductance G r ≈ α r A 0 .In this limiting case the heat transfer between the bodies is determined by the propagating electromagnetic (EM) waves radiated by the bodies and does not depend on the separation d between the bodies.Electromagnetic waves (or photons) always exist outside any body due to thermal or quantum fluctuations of the current density inside the body.The EM-field created by the fluctuating current density exists also in the form of evanescent waves, which are damped exponentially with the distance away from the surface of the body.For an isolated body, the evanescent waves do not give a contribution to the energy radiation.However, for two solids separated by d < d T , the heat transfer may increase by many orders of magnitude due to the evanescent EM-waves; this is often referred to as photon tunneling.
For short separation between two solids with flat surfaces (d << d T ) the heat current due to the evanescent EM-waves is given by [15,18,19] where where ǫ(ω) is the dielectric function, orientations average if anisotropic.We assume again T 1 = T 0 +∆T with ∆T << T 0 .Expanding (6) to linear order in ∆T gives J = α e ∆T with where η = ̵ hω k B T and ξ = qd.Note that the second integral depends on ω so the first integral involves also the result of the second integral.
From (7) it follows that the heat current scale as 1 d 2 with the separation between the solid surfaces and we write α e = a(T ) d 2 .For a spherical particle in contact with another particle (or a flat surface) the interfacial separation d ≈ r 2 2R (where 1 R = 1 R 0 + 1 R 1 ) where r is the radial distance in the xy-plane from the center of the contact region.In this case it is trivial to obtain the heat conductance where d 0 ≈ 0.3 nm is the equilibrium separation between the surfaces at the contact point.Consider now two clean surfaces of (amorphous) silicon dioxide (SiO 2 ).The optical properties of this material can be described using an oscillator model [20]  The a(T ) factor in the heat transfer coefficient αe = a d 2 as a function of temperature for two silica surfaces (separation d).The red line is obtained using the two oscillator model (8) while the green line is obtained using the numerical data for the dielectric function given in Ref. [20].
The frequency dependent term in this expression is due to optical phonons.The values for the parameters ǫ ∞ , (a 1 , ω 1 , γ 1 ) and (a 2 , ω 2 , γ 2 ) are given in Ref. [20].In Fig. 7 we show the temperature dependency of the a(T ) parameter obtained using (8) (red line) and also as obtained using directly the measured ǫ(ω) (green line).For T = 200 K we get a(T ) ≈ 1.0 × 10 −12 W K and we use this value below when calculating the contact conductance G.In Appendix C we show the a(T ) function for another mineral (olivine) of interest in applications to asteroids.
In the present case the heat transfer is associated with thermally excited optical (surface) phonons.That is, the electric field of a thermally excited optical phonon in one solid excites an optical phonon in the other solid, leading to energy transfer.The excitation transfer occur in both u(x,y) rigid solid rigid solid J(x,y)~1/u 2 FIG. 8.The heat current between the two surfaces is assumed to depend on the separation as J(x, y) ∼ u −2 (x, y).directions but if one solid is hotter than the other, there will be a net transfer of energy from the hotter to the colder solid.

Radiative contribution to G: numerical results
We have calculated the heat contact conductance using the same approach as used to study the adhesion between particles with random surface roughness in Ref. [11].No two natural stone particle have the same surface roughness, and the heat transfer between two particles will depend on the particles used.To take this into account we have generated particles (with linear size L = 2R) with different random surface roughness but with the same surface roughness power spectrum.That is, we use different realizations of the particle surface roughness but with the same statistical properties.For each particle size we have generated 60 particles using different set of random numbers.The surface roughness was generated as described in Appendix A of Ref. [21] by adding plane waves with random phases φ q and with the amplitudes determined by the power spectrum: q⋅x+φq)   where B q = (2π L)[C(q)] 1 2 .We assume isotropic roughness so B q and C(q) only depend on the magnitude of the wavevector q.
We have used nominally spherical particles with 6 different radii, where the radius increasing in steps of a factor of 2 from R = 78 nm to R = 2.53 µm.The longest wavelength roughness which can occur on a particle with radius R is λ ≈ 2R so when producing the roughness on a particle we only include the part of the power spectrum between q 0 < q < q 1 where q 0 = π R and where q 1 is a short distance cut-off corresponding to atomic dimension (we use q 1 = 1.4 × 10 10 m −1 ).We will refer to these particles as granite particles because the power spectra used are linear extrapolation to larger wavenumber of the measured granite power spectrum.For more details about the numerical procedure see Ref. [11].
We will now present numerical results for the heat conductance of granite particles.We will also consider particles with the same sizes as above but with larger and smaller surface roughness, obtained by scaling the height h(x, y) for the granite particles with scaling factors s = 0 (smooth surface), 0.1, 0.3, 1 and 2. Note that scaling h(x, y) by a factor of s will scale the power spectrum with a factor of s 2 but it will not change the slope of the C(q) relation on the log-log scale so the Hurst exponent (and the fractal dimension) is unchanged.
We assume that the heat current depends on the surface separation u(x, y) as J(x, y) ∼ u −2 (x, y).This holds accurately only in the small slope approximation.The heat conductance G e is obtained by integration α e = a(T ) u 2 (x, y) over the surface area (see Fig. 8).
Fig. 9 shows the cumulative probability for the conductance G for all the particles with different radius.The probability distributions are obtained by using for each particle size 60 different surface roughness realizations with the same power spectra.The calculations are for the granite surface (scaling factor s = 1).Note that there is a slight increase in G with increasing particle radius.For the Van der Waals interaction between the particles the different particle radius gave nearly the same cumulative probability distribution i.e., the pull-off force, and the statistical fluctuations in the pull-off force, where nearly the same for all the particles.The reason for why G exhibit a stronger (but still very weak) dependency on the particle radius is that while G depends (for flat parallel surfaces) on the interfacial separation d as 1 d 2 while the Van der Waals interaction decay faster as 1 d 3 ; see, e.g., [42].
Fig. 10 shows the heat contact conductance as a function of the particle radius (log-log scale) at T = 273 K.The solid lines are assuming the evanescent-wave electromagnetic (EM) coupling between the particles.The The heat contact conductance as a function of the particle radius (log-log scale) at T = 200 K, H=1, parameter s.The solid lines are assuming the evanescent-wave electromagnetic coupling between the particles.The dashed line is the result assuming propagating-wave electromagnetic coupling (i.e.black body radiation) but this result is only for illustrative purpose because the Stefan-Boltzmann law is not valid for the small surface separation occurring with particles as small as studied here.The different solid lines are for particles where the surface roughness of the granite particle is scaled with different factors s between 0 (smooth surface) and 2. The heat contact conductance as a function of the particle radius (log-log scale) at T = 200 K for particles with self-affine fractal surface roughness with different Hurst exponent H = 0.6, 0.7, 0.8 and 0.9; roughness amplitude s=1.All surfaces have the same root-mean-square slope when including the roughness on all length scales (see Ref. [11] for more details about the surface roughness).
dashed line is the result assuming propagating EM-wave coupling (i.e.black body radiation) but this result is only for illustrative purpose because the Stefan-Boltzmann law is not valid for the small surface separation occurring with particles as small as studied here.The different solid lines are for particles where the surface roughness of the The cumulative heat conductance obtained by integration the heat current ∂J(x, y, T ) ∂T over a circular region x − x0 < r with the radius r, and centered at the point x0 where the heat current is maximal.For one realization of the surface roughness for all the particles and T = 273 ○ C.
granite particle is scaled with different factors between 0 (smooth surface) and 2. Note that for particles with radius R > 10 µm the black-body radiation will dominate the heat transfer for granite particles but for particles with smooth surfaces (s = 0) a much larger particle radius is needed before the black body radiation dominates the heat transfer.Note also that the surface separation d must obey d > d T = c ̵ h k B T ≈ 10 µm (for T = 273 K) in order for the Stefan-Boltzmann law to be valid.Fig. 11 shows the heat conductance as a function of the particle radius (log-log scale) at T = 273 K for particles with self-affine fractal surface roughness with different Hurst exponent H = 0.6, 0.7, 0.8 and 0.9.All surfaces have the same root-mean-square slope when including the roughness on all length scales (see Ref. [11] for more details about the surface roughness).
Fig. 12 shows the spatial dependency of the heat current ∂J(x, y, T ) ∂T close to the point where it is maximal.The results is for one realization of the surface roughness for the particle with the radius R = 0.32 µm and T = 273 ○ C. From the figure it may appear that most heat flow is localized to a few nanometer sized region close to the point where the current in maximal.However, when plotted on a logarithmic scale it is clear that the biggest contribution to the heat transfer occur from a much bigger surface region.Thus, in Fig. 13 we show the heat current ∂J(x, y, T ) ∂T along the x and y-directions through the point where it is maximal for a particle with the radius R = 2.53 µm.In Fig. 14 we show for all the particles the cumulative heat conductance obtained by integration the heat current ∂J(x, y, T ) ∂T over a circular region x − x 0 < r with the radius r, and centered at the point x 0 where the heat current is maximal.Note that for the 3 biggest particles about half of the total heat conductance result from the heat flow within a circular region with the radius r ≈ 0.1 µm.This is an important result because it implies that the constriction contribution to the thermal resistance for the EM heat transfer can be neglected.Thus if we assume r e ≈ 0.1 µm we get the constriction resistance 1 (2K 0 r e ) ≈ 10 7 K W which is much smaller than the resistance R e = 1 G e ≈ 10 10 K W resulting from the evanescent EM waves.Since these resistance act in series we can neglect the constriction contribution.Hence for very rough particles as produced by crunching mineral stone one expect the contact conductance to be determined by the evanescent EM waves.This conclusion is summarized in Fig. 15.
In the calculation above we have assumed that both materials in the contact region between two particles are identical.Many mineral particles consist of grains with (slightly) different chemical composition.In that case the two reflection factors R 0 (ω) and R 1 (ω) in (7) will differ, which will reduce a(T ) as the heat transfer depends on the product ImR 0 (ω)ImR 1 (ω), which depends on Imǫ 0 (ω)Imǫ 1 (ω).
For completeness me mention the "grain" thermal con- ductance, G p , which can be estimated by approximating the particle (material thermal conductivity K 0 (T )) by a cube with equal volume, leading to G p in W/K is of the order of the grain diameter in m, in practice, 10 −6 to 10 −2 W/K, so small that it can in almost all cases be neglected.Note that the classical constriction resistance is based on continuum theory and neglects the phonon interface (acoustic mismatch) resistance 1 G I F , which exists even for strong coupling if the bonded lattices are even slightly different or have a different orientation (usually this is the case at grain boundaries).A typical value [28] is 10 8 W m 2 K for the conductance per area, so e.g. for 10 µm grains, G I F ≈ 0.01 W K.

Comparison with experiment: a proposal
The theory presented above could be tested experimentally by measuring the heat conductivity of granular materials under high vacuum condition.We suggest to use crunched minerals sieved to obtain monodisperse particle size distribution.We suggest to use crunched pure minerals (crystalline preferred), sieved to obtain a monodisperse particle size distribution.The roughness power spectrum of the mineral particles shall be determined and their dielectric function in the mid-IR must be known.Thermal conductivity measurements should be done in vacuum, after baking out the samples to make sure the particles lost all their physisorbed water.Caution ahould be exercised w.r.t.boundary conditions in all experiments, whether it be a line heat source, guarded hot plate, strip heat source, or whatever: the packing structure (and contact resistance) is always disturbed near a boundary or implanted wire/strip.Measurements should be performed for particle sizes from ∼ 1 µm to at least ∼ 500 µm such that the transition from dominating near-field EM heat transfer to, additionally, firmly geometrics optics radiative far-field conductivity can be observed.The temperature range should be as large as possible, to confirm the temperature dependence of a(T ) and to be able to separate the contact conductance from the radiative (far field) conductivity, which ought to be possible for very low (say 80 K) and very high (say 600 K) temperatures, where dlog(a) dlogT typically has slopes quite different from 3. Temperatures down to 25 K with e.g.crystalline corundum (synthetic sapphire) Al 2 O 3 would be interesting too, since this crystal has a very large (200 W/K/m) bulk thermal conductivity peak at 30 K, two orders of magnitude higher than at room temperature [41]; this should be clearly seen in the data if solid conduction were phonon-dominanted via strong contacts.Various porosities from random loose to random close packing should be prepared by a dedicated tapping procedure and determined precisely.It would also be interesting to perform the same experiments in the normal atmosphere at different humidities to include capillary bridges and heat diffusion in the fluid and in the surrounding air [6].Additionally, lithostatic pressure could be applied to study the its effect on contact conduction.To analyze the latter experiments one would need to extend the study we present in this paper.
Interesting heat conductivity data was presented by Sakatani et al. [8].They used spherical glass beads with radii ranging from 2.5 µm to 427 µm.The roughness of the glass beads was not studied in detail but SEM images showed non-random roughness.The experiments was performed in vacuum, but liquid bridges may have formed while handling the granular material in humid lab air before the vacuum study.Experiments have shown that under the influence of water or humidity Si-O-Si bonds may form between two silica surfaces in the contact regions which could remain in the vacuum condition and reduce the interfacial contact resistance [22][23][24][25].However, if the granular media is exposed to vibrations when in the vacuum then these chemical bonds may be broken in which case the particle-particle interaction may be mainly of the Van der Waals type as assumed here.For these reason no quantitative comparison of the data of Sakatani et al. with the present theory is possible or meaningful.

Discussion
Most asteroids larger than ∼ 100 m are rubble piles, i.e. consist of more than one solid objects, and are usually covered by a regolith layer (dust, broken rocks, and other related materials).In fact all asteroids which have been studied with spacecraft consist of a wide distribution of fragments with sizes ranging from about 100 meters to centimeter or less.The average temperature in the asteroids (typically ≈ 200 K) is so high that on the time scale of billions of years all loosely bound molecules have desorbed so one expect no mobile adsorbed molecules which could form capillary bridges between the particles (or fragments).Since the gravitational field is too weak to allow for an atmosphere the particles in asteroids are surrounded by vacuum.Hence the heat transfer between the particles can only occur via the area of real contact The a(T ) factor as a function of temperature (from Fig. 7) on a log-log scale (solid line).The dotted line has a slope of 3 corresponding to a temperature dependency T 3 expected from the black-body heat transfer.
For very open (low density) granular material photons can travel long distances resulting in high effective thermal conductivity photon T 1 T 2 FIG. 17.In very open (low filling factor) granular materials cavities much bigger than the particle diameter can occur.In this case (Stefan-Boltzmann) heat radiation may result in large effective thermal conductivity.
(which are nanosized regions as shown in Ref. [11]) or via the fluctuating EM-field (propagating or evanescent) which occur around all solid objects.We have shown that because of the weak adhesion and gravitational force, the contribution to the contact conductance from the area of real contact is negligible compared to the contribution from the EM-field.
From studies of the surface temperature of asteroids one can deduce the thermal inertia Γ = (ρCK) 1 2 .Most asteroids have thermal inertia between 10 and 1000 Jm −2 K −1 s −1 2 and typically Γ ≈ 100 Jm −2 K −1 s −1 2 (see Ref. [10]).Assuming ρ ≈ 2 × 10 3 kg m 3 , C ≈ 500 J kgK and Γ ≈ 100 Jm −2 K −1 s −1 2 we get K ≈ 10 −2 W Km. This is much smaller than the thermal conductivity of silicate rocks which is ≈ 1 − 10 W Km. This proves that most asteroids does not consist of single silica blocks but of weakly coupled fragments where most of the thermal resistance comes from the thermal coupling between the fragments.
Assume that all the particles in an asteroid would have equal size.In this case it is interesting to consider two limiting cases where the heat transfer occur either via the radiative (black-body) EM-field, or via the evanescent EM-field.In the second case G is nearly independent of the particle radius and of order G ≈ 10 −9 W K giving R ≈ G K ≈ 0.1 µm.In the first case G = µR 2 where µ ≈ 5.7 W m 2 K for T = 200 K.This gives R = K µ ≈ 2 mm.In reality there is a wide distribution of particle sizes and the heat transfer may involve both the radiative and the evanescent EM-wave interaction.However, if the big particles (fragments) are surrounded by a matrix of much smaller (micrometer-sized) particles (dust), as suggested by other studies [11], then the radiative coupling my be less important than the evanescent contribution.In particular, unless the cavity regions between the particles are of size d > d T = c ̵ h k B T ≈ 10 µm the the Stefan-Boltzmann law will not be valid.In this case the big particles will effectively act as regions of infinite heat conductivity (because heat transfer in compact solids fast), and the resistance to the heat flow will occur mainly in the (low volume) matrix of small particles surrounding the bigger particles.Thus taking this effect into account, and assuming that the small particles have radius R ≈ 1 − 10 µm as suggested by other studies [11], result in thermal inertia similar to what is observed.Fig. 16 shows the a(T ) factor as a function of temperature (from Fig. 7) on a log-log scale (solid line).The dotted line has a slope of 3 corresponding to the temperature dependency T 3 expected from the black-body heat transfer.Fig. 16 shows that the contribution to the heat transfer from the evanescent waves may increase faster or slower with temperature than the black body contribution, depending on the temperature region.
We note that in very open (high porosity) granular materials cavities much bigger than the particle diameter can occur (see Fig. 17).In this case (Stefan-Boltzmann) heat radiation may result in large effective thermal conductivity.The (macro-)porosity φ [which equals (void volume)/(total volume)] of rubble pile asteroids ranges from φ ≈ 0.15 (see Ref. [12]) to φ ≈ 0.5 (see Ref. [26]), and the very top layer of lunar regolith has a porosity > 0.8 (see [27]) which may result in some fraction of the emitted photons traveling one or several particle diameters.
In the asteroid research field it is usually assumed that the constriction resistance is a very important contribution to the thermal contact resistance.We have shown that the constriction resistance is important only when cold-welded (for metal) or sintered contact forms, or where strong (chemical) bonds form between the particles in the contact area.However, experiments have shown that the asteroid particles interact only with very weak forces and the contact resistance is dominated by the interfacial thermal resistance (determined by ( 2)), and the constriction resistance can be neglected.
In asteroid research it is usually assumed that the thermal conductivity K eff = A + BT 3 where A is assumed to be due to the constriction resistance and hence proportional to the bulk thermal conductivity K 0 (T ) which is only weakly temperature dependent.However, we find that the constriction resistance term can be neglected and that A may be due to the evanescent EMwaves K e ∝ a(T ), where a(T ) typically has a strong T -dependence, where at intermediate temperatures it is roughly ∝ T 3 .Thus, K e may have a similar Tdependence as K r .The particle size-dependence of K e is also fundamentally different, K e ∝ (1 R)R n , where n ≈ 0.2, while in the conventional theories the term A is independent of R, or increases only for very small R due to an (assumed) larger effect of van Der Waals adhesion deformation.This should be observable in data for small R <≈ 10µm.
Finally we note that the Hertz-or JKR-theory for the contact radii between two particles in the conventional theories is unrealistic.The particles have large roughness and are made from stiff materials and the JKR (and Hertz) theory is not valid.Instead the adhesive interaction is very weak as we already have shown [11].This also means that the theoretical dependence of K on the external (e.g., lithostatic) pressure P , roughly ∝ P 1 3 (plus a constant for P = 0), is not due to the contact mechanics between smooth adhering spherical particles.What matters, is the reduction of the surface separations of the two rough particles in their contact zone upon an externally applied pressure, since this controls evanescent EM-wave heat transfer.

Acknowledgments
Part of this work was performed when Bo Persson was participating in the UCSB, KITP program "Emerging Regimes and Implications of Quantum and Thermal Fluctuational Electrodynamics".This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.
Appendix A In Sec. 3 we considered heat diffusion in a simple cubic arrangement of spherical particles with radius R and found K = γG R with γ = 3 π.The study assumed that the heat transfer occur only at, or in the vicinity of, the contact regions between the particles as expected for the phononic contribution K a or the contribution from the evanescent EM waves K e .In Sec. 4 we also considered the radiative (black-body) contribution to the heat conductance G r ≈ πR 2 σT 3 so we expect K r = βσRT 3 where β ≈ 3.Here we will give approximate expressions for γ and β for realistic distribution of particles as derived by Arakawa et al [29] and by Ryan et al [30], respectively.We expect β and γ to be of order unity, but they depend on the particle filling factor (or porosity), and reliable estimates of these factors are needed for comparison with The geometrical factor γ for heat diffusion via particle contacts, as a function of porosity (for random packings).After [29] (solid line) and [8] but with the coordination number model of [29] (dashed line).The value for the regular cubic packing is indicated, too.
experiments or observational data.

Contact conductivity K
For a general granular medium we have to modify equation ( 1) by a factor which depends on the porosity, and on packing geometry (ordered or random; shape and friction of particles).
Arakawa et al [29] have performed simulations of the heat transfer in granular media consisting of micrometersized spherical particles with random packing with the porosity ranging from 0.37 to 0.99.The simulated results where fitted to obtain the effective thermal conductivity K = γG R with where the average coordination number C depends on the porosity φ.For spherical particles 62 , but for irregular particles C is typically smaller than for spheres at the same porosity [31,32].For a simple cubic lattice φ = 1 − π 6 ≈ 0.476 giving C ≈ 5.3 and γ ≈ 0.31 which is smaller than the factor 3 π ≈ 0.95 expected for a simple cubic lattice of spherical particles.The factor γ is shown as a function of the porosity in Fig. 18 where the square symbol indicate γ for the simple cubic lattice.
Radiative conductivity K r Ryan et al. [30] have shown that for opaque particles in the geometrical optics regime K r = βσRT 3 with FIG.19.The geometrical factor β for conventional radiative conductivity Kr = βσRT 3 , as a function of porosity for emissivity ǫ = 1 in the limit Λ → ∞.Adapted from [30].
For small enough spherical particles, or particles with high enough thermal conductivity K 0 , the temperature in the particles will be uniform corresponding to Λ → ∞.In this limit f ≈ 1.For a simple cubic lattice of particles φ ≈ 0.476 and with ǫ = 1 we get F ≈ 1.15 giving β ≈ 9.2.Fig. 19 shows β as a function of φ when Λ = ∞.
Equation (A2) is strictly valid only in the regime of geometrical optics, i.e., if the thermal wavelength (in meter) 2.9 × 10 −3 T −1 (with T in Kelvin) is larger than particle diameter D. The case when the thermal wavelength is smaller than the particle size D and/or the particles cannot be regarded as opaque to thermal radiation was studied in Ref. [34,[37][38][39].However, in this case the contribution from the evanescent EM waves may be more important.
To facilitate the application of our new heat transfer model for the planetary sciences or quantitive comparisons with experimental data, we need correlation equations that represent G e as a function of R, s, H and then the temperature-dependence (T 0 = 200K): As for G e (R, s, H) silica−silica, T0 it is possible to interpolate the results shown in figures 10, 11.Of course, the applicable values for Hurst exponent H and roughness scale factor s have to be known or estimated.For comparison, we state these properties for crushed meteorite fragments from previous work [11, fig. 11] using the data by Nagaashi et al. [43] of meteorite fragments of type CM2, CV3, LL3.5 and glass beads.We found s = 1.36 ± 0.11 for Murchison meteorite (CM2) fragments, s = 3.0 ± 0.2 for Allende meteorite (CV3) fragments, s = 4.3 ± 0.4 for NWA 539 (LL3.5)fragments, and, tentatively, s = 0.23 ± 0.03 for the glass beads (GB); this assumes that H = 1 (exact) which is consistent with the data and common for natural surfaces.
The temperature dependence of G e is only in a(T ) which depends significantly on the combination of minerals across the contact, i.e., their dielectric functions in equation (7).Dislike minerals generally produce smaller a.For the major rock-forming minerals and other compounds of interest, i.e., pyroxenes, feldspars, phyllosilicates, olivines, ice Ih, quartz, calcite, dolomite, magnesite, silica and industrial glasses, dielectric functions are readily available in the literature (usually given as complex refractive index for the near-and mid-IR) and the integral ( 7) is easy to evaluate.In a typical rock (fragment), the constituting minerals (with volume fractions assumed to be known) are spatially distributed as grains, lamelleae etc. on many length-scales.It will be the subject of future work to find a proper weighted average of the a(T )-factors of every possible combination.

Appendix B
Let F (t) be the pulsating force exerted on the solid 1 from an asperity on solid 0. We assume the force localized to the point x = 0 on the surface of solid 1.Thus the stress act on the solid 1.We get The energy transfer to solid 1 from the pulsating force Using that [35] u z (q, ω) = M zz (q, ω)σ z (q, ω) Next if we write q = (ω c T )ξ (see Ref. [35]) then Since ∆E is real we can write Thus we can write We assume the fluctuation in the position of the atoms are small so that we can expand the interaction potential U (d) at the asperity contact in the displacement away from the equilibrium separation d = d eq to linear order in the displacement s = d − d eq so that where k = U ′′ (d eq ).Thus we get Ḟ (t) = kv z (t) where v z = ṡ(t) is the velocity of the atom normal to the surface.Thus from (B2) we get We can write dt v 2 z (t) = t 0 ⟨v 2 z ⟩ where t 0 is the time period for which we calculate the energy transfer.Thus the energy transfer per unit time where ⟨v 2 z ⟩ is the time average (or ensemble average) of the fluctuating atom velocity v z (t).We first assume high temperatures so that we get G a = k B η.We have shown in Ref. [13,35] that The derivation above assumes high temperatures, but a result valid for all temperatures can be obtained as follows: Note that phonons are harmonic oscillators where on the average half of the total energy is potential (elastic deformation) energy and half is kinetic energy.Thus the total (phononic) energy E in a solid with N atoms is Thus the heat transfer will depend on where C V is the heat capacity.Thus the heat conductance G a = k B η is valid for all temperatures if we replace k B with C V 3N .For high temperatures C V ≈ 3N k B and we recover the limit considered above.The physical picture behind the heat transfer described above is that the irregular (random) thermal movement of the atoms in the contact region (here of atomic size) exert pulsating forces on the opposite solid (or particle) which result in phonon emission.This emission of phonons occurs in both direction but is stronger from the hotter to the colder solid.It is assumed the interaction between the two solids in the contact region (as manifested by the spring constant k) is so weak that it does not influence the irregular motion of the atoms.Here it is interesting to note that the heat capacity can be expressed as an energy fluctuation term so the heat transfer conductance is the product of a friction coefficient (or inverse relaxation time) and an energy fluctuation term.Using the Debye model, where the phonon dispersion in the bulk is assumed linear for all phonon wavenumbers up to a cut-off wavenumber ω D , one gets where T D = ̵ hω D k B and x D = T D T .The heat capacity for silica is well described by the Debye model with T D ≈ 364 K (see Ref. [36]).For T = 200 K we get T T D ≈ 0.55 and for this relative temperature the heat capacity has already reached 85% of it high temperature value.Hence using the high temperature expression (2) for the damping rate is a good approximation.which depends on Imǫ 0 (ω)Imǫ 1 (ω).
Fig. 20 shows the imaginary part of the dielectric function of silica (red line) and olivine (blue) as a function of the frequency in the infrared region.In all cases the loss function Imǫ(ω) is due to optical phonons.Fig. 21 shows the a(T ) factor in the heat transfer coefficient α e = a d 2 as a function of temperature for two flat surfaces (separation d) of silica (red line) and olivine (blue), and for a silica and olivin surface (green).

FIG. 1 .
FIG. 1.(a) Optical image of the asteroid Bennu.(b) The temperature distribution obtained from Bennu thermal emission.Credit: NASA/Goddard/University of Arizona

FIG. 2 .
FIG.2.Heat flow in a simple cubic lattice of spherical particles (radius R).We assume the temperature depends only on one coordinate direction denoted by x.Random packings of various porositys are considered in appendix C.

FIG. 6 .
FIG. 6.The imaginary part of the dielectric function of silica in the infrared region.The two high peaks are due to two different optical phonons.

FIG. 9 .
FIG.9.The cumulative probability for the evanescent EM contribution to the contact conductance G.The probability distributions are obtained from 60 simulations for each particle radius.The 60 simulations use 60 different realizations of the particle surfaces topography but with the same power spectra.The calculations are for the granite surface (scaling factor s=1, Hurst exponent H=1) at 200 K.
FIG. 10.The heat contact conductance as a function of the particle radius (log-log scale) at T = 200 K, H=1, parameter s.The solid lines are assuming the evanescent-wave electromagnetic coupling between the particles.The dashed line is the result assuming propagating-wave electromagnetic coupling (i.e.black body radiation) but this result is only for illustrative purpose because the Stefan-Boltzmann law is not valid for the small surface separation occurring with particles as small as studied here.The different solid lines are for particles where the surface roughness of the granite particle is scaled with different factors s between 0 (smooth surface) and 2.
FIG. 11.The heat contact conductance as a function of the particle radius (log-log scale) at T = 200 K for particles with self-affine fractal surface roughness with different Hurst exponent H = 0.6, 0.7, 0.8 and 0.9; roughness amplitude s=1.All surfaces have the same root-mean-square slope when including the roughness on all length scales (see Ref.[11] for more details about the surface roughness).

FIG. 13 .
FIG.12.The heat current ∂J(x, y, T ) ∂T close to the point where it is maximal.For one realization of the surface roughness for R = 0.32 µm and T = 273 ○ C.

FIG. 15 .
FIG.15.The EM (evanicent waves) contact resistance and the constriction resistance act in series and since the former is much larger the latter can be neglected.
FIG. 18.The geometrical factor γ for heat diffusion via particle contacts, as a function of porosity (for random packings).After[29] (solid line) and[8] but with the coordination number model of[29] (dashed line).The value for the regular cubic packing is indicated, too.

A
similar expression with T 0 replaced by the temperature T 1 gives the energy transfer Q1 from solid 1 to solid 0. Thus the net energy transfer isQ = Q0 − Q1 = G a (T 0 − T 1 )

FIG. 20 .FIG. 21 .
FIG.20.The imaginary part of the dielectric function of silica (red line) and olivine (blue) as a function of frequency in the infrared region.The loss function structures are due to optical phonons.Data from Ref.[20] and[40].