Subscriber Authentication Point
Free Access
 Issue A&A Volume 519, September 2010 A75 16 Planets and planetary systems https://doi.org/10.1051/0004-6361/201014281 15 September 2010
A&A 519, A75 (2010)

## Thermal stresses in small meteoroids

D. Capek1 - D. Vokrouhlický2

1 - Astronomical Institute of the Academy of Sciences, Fricova 298, 251 65 Ondrejov, Czech Republic
2 - Institute of Astronomy, Charles University, V Holesovickách 2, 180 00 Prague 8, Czech Republic

Received 18 February 2010 / Accepted 24 April 2010

Abstract
Aims. We evaluate thermal stresses in small, spherical, and homogeneous meteoroids with elastic rheology and regular rotation. The temperature variations are caused by the absorbed sunlight energy being conducted into the interior layers of the body. Our model assumes arbitrary thermal conductivity value, but restricts itself to a linearized treatment of the boundary conditions of the heat diffusion problem. We consider the diurnal insolation cycle only as if the body were in a fixed position along its heliocentric orbit. This constrains the upper limit to the object size to which our modeling is applicable.
Methods. We derive analytical expressions for the components of the thermal stress tensor throughout the body. Using two sets of material properties (ordinary and carbonaceous chondrites), we study the conditions required for material failure caused by thermal stress leading to fission.
Results. Our results indicate that the onset of thermal failure in the meteoroid depends on a number of parameters including the heliocentric distance, the size, the rotation frequency, and the orientation of the spin axis with respect to the solar direction. In our case, we find large, centimeter- to meter-size, slowly rotating meteoroids or those with a spin axis pointing towards the Sun or both, are the most susceptible to the thermal bursting. This may have implications for the (i) size distribution of meteoroids in various streams depending on their heliocentric orbit and the physical characteristics of their parent bodies; (ii) orbital distribution of sporadic complexes of meteoroids in the planet-crossing zone; and/or (iii) fate of fragments released during comet disintegration events, especially those with low perihelia (e.g., Kreutz class).

Key words: meteorites, meteors, meteoroids - minor planets, asteroids: general - methods: analytical

## 1 Introduction

Small bodies in the inner parts of the Solar System are frequently brought to small heliocentric distances, well within the orbit of Mercury. They might either originate in orbits with such a low pericenter value (e.g., particles in meteoroid streams; Jenniskens 2006) or temporarily reside in these orbits during an evolution driven by chaotic dynamics in the planet-crossing region (e.g., meteoroids or small near-Earth asteroids; Marchi et al. 2009). In the harsh radiation environment within a distance of 0.15 AU from the Sun (equilibrium temperatures 700 K), the surface layer and/or the whole volume, for objects of small sizes, can undergo interesting physical alteration processes such as release of the volatile elements and/or partial melting of silicate components (e.g., Marchi et al. 2009; Capek & Borovicka 2009). Those effects can directly influence spectral characteristics, put cometary objects into either a depleted or dormant state, and/or accelerate space weathering processes.

Yet another property related to the strong solar heating at small heliocentric distances is the growth of the thermal stress, namely a mechanical stress caused by a non-uniform temperature field in the solid body that, under some circumstances, may exceed the material strength. Either crack formation or a break-up of the body may occur (thermal shock material failure). Mechanical integrity of the material may be affected even when the thermal stress value does not directly reach the critical failure value. This is the case for periodic temperature variations with large enough gradients that produce a slow but steady propagation of pre-existing micro-cracks leading eventually to the break-up. In this case, we refer to the thermal fatigue (Hall 1999).

The role of thermal stresses in the physics of small Solar System bodies has been studied by many of authors. Their importance to the evolution of cometary nuclei was advocated and discussed, for example, by Kuehrt (1984), Tauber & Kuhrt (1987), and Tambovtseva & Shestakova (1999). In their view, as the cometary nucleus approaches the Sun, the thermal stresses may become several orders of magnitude larger then the estimated solar tidal stresses. They may be able to exceed the material strength locally or globally, producing cracks at the surface, resulting in an increase in the cometary activity, or leading to a splitting of the whole nucleus.

Grinin et al. (1996) pointed out that even non-icy planetesimals residing on highly eccentric orbits can be destroyed by thermal stresses in the neighborhood of UX Orionis-type young stars producing a heavy-element-rich gaseous envelope radiatively expelled from the parent region. Shestakova & Tambovtseva (1997) applied similar ideas to our own Solar System. In their scenario, meteoroids approaching the Sun on highly eccentric orbits experience gradual heating. A steep thermal gradient between the cold core and hot surface slowly grows up, resulting in thermal stresses that can eventually result in thermal fission for bodies of range from about 10 cm to 10 m. For larger bodies, in their model, the thermal gradient and the associated thermal stresses are not large enough to cause material failure. These authors also noted that surface melting at very high temperatures may be able to stop the steepening of the thermal gradient between the surface and the core. At lower temperatures, sublimation may indeed lead to surface cooling. Another important process may be the built-up of a low-conductivity surface (regolith-like) layer that may shield the bulk of the body from large temperature variations and thus thermal stresses. This is most likely the fate of large bodies that may be able to form this protective layer on their surface. On the other hand, small bodies are heated throughout their volume by efficient thermal conduction and become nearly isothermal. The lack, or small value, of the thermal gradient in this case means that the associated thermal stresses are insufficient for break-up to occur.

While the previous analyses of Kuehrt (1984), Shestakova & Tambovtseva (1997), and Tambovtseva & Shestakova (1999) considered the time dependence of the temperature field caused by the approach to the Sun along a parabolic orbit, they adopted significant simplifications that our approach attempts to remove. Most importantly, they (i) circumvented the difficulty in solving the thermal boundary condition by assuming exact equilibrium between the incident solar flux and thermal emission (neglecting heat conduction into the body); and (ii) assumed only a radial temperature dependence (neglecting the latitudinal and longitudinal dependence related to the specific location of the Sun). Thus, according to their approximation, the thermal stresses result merely from temperature difference between the cold-maintained core and the solar-heated surface. But in reality the temperature is also affected by a variation in the insolation caused by the body's rotation about the spin axis. While bodies with rapid rotation typically have axially symmetric temperature profiles with only latitudinal temperature gradients complementing the radial ones, those that can rotate slowly can also build longitudinal thermal fields and gradients. In that case, the overall analysis of the thermal field may be complicated, but may hold important information about whether the conditions of thermal failure originate near the center of the body or affect subsurface layers. With this goal in mind, we refrain from including seasonal temperature variations in this work, leaving this to a forthcoming paper, and focus on including the whole complexity of non-radial features in the temperature field. This is the main novelty of this work compared to previous studies mentioned above.

As discussed in some detail in the next section, the lack of seasonal features in the time dependence directly causes a size constraint to which our approach is applicable. It turns out, though, that our results are fully applicable to typical particles in the meteoroid streams (sizes 10 cm, say), which are the main focus of this paper. By using analytical, rather than numerical, methods, we trade quantitative exactness for qualitative understanding. This obviously requires necessity to adopt simplifications, most of which are summarized in the next section.

## 2 Theory

### 2.1 Formulation of the problem

Our goal is to determine thermal stress field in small meteoroids using an analytical approach and analyze the conditions required for their thermal disruption. To ensure, that our calculations are manageable, we adopt a number of simplifying assumptions. The generally irregular shape of meteoroids is described by a sphere of radius R. In the same way, the generally tumbling state of rotation is represented by regular rotation about a fixed (spin axis) direction and a constant rotation angular frequency . Based on these assumptions, it is most straightforward to describe physical quantities, such as the temperature or stress tensor fields, using spherical coordinates with the origin r=0 at the center of the body, colatitude measured from the spin axis , and longitude being defined arbitrarily in our approach. Physical parameters, such as thermal and elastic constants, are assumed to be homogeneous, isotropic, and independent of the temperature. We also restrict ourselves to the approximation of elastic rheology.

### 2.2 Temperature distribution

As far as the temperature T distribution in the body is concerned, we assume that it responds primarily to the solar heating. This means that we assume quasi-static thermoelasticity (Parkus 1976), neglecting the deformation field as a heat source. In addition, we adopt two assumptions about the meteoroid rotation frequency and size: (i) the former is low enough for the thermal relaxation to occur throughout its whole volume on a much shorter timescale than the revolution period about the Sun. For reasonable material parameters, the thermal relaxation timescale is not much longer that the rotation period. Therefore, the instantaneous heliocentric position determines the thermal and stress state of the meteoroid independently of its orbital history; (ii) the diurnal-cycle approximation from (i) has also a consequence on the assumed maximum size of the meteoroid. This is because the modeled insolation field has all mean-motion-induced Fourier modes mistakenly'' collapsed to zero frequency. This static part of the temperature field can penetrate arbitrarily deeply into the body, while the true minimum frequency of the insolation energy loading, namely the mean motion frequency n of the heliocentric revolution, allows the temperature field to penetrate only to a certain maximum depth (usually called the penetration depth of the seasonal thermal wave). One can easily see that (e.g., Vokrouhlický et al. 2007; Vokrouhlický 1999)

 (1)

where Kis the thermal conductivity, is the bulk density, and cis the heat capacity. Our results in this paper are thus valid only for bodies with . For characteristic values of meteoroids (see e.g. Sect. 2.5), we obtain R < 1-4 m. This constraint is to be kept in mind especially when discussing the implications of our results in Sect. 4.

Having defined the approach and simplifying assumptions, we now proceed with the solution for the temperature distribution. This is determined by solving the heat diffusion equation

 (2)

complemented with appropriate boundary conditions, which in our case are (i) regularity of the solution in the center; and (ii) energy conservation at the surface. The latter can be expressed as

 (3)

where is the thermal emissivity,  W m-2 K-4 is the Stefan-Boltzmann constant, is the unit outer normal to the surface, A represents the albedo value, and is the solar radiation flux. In principle, A is the hemispheric albedo value, which depends on the zenith angle of the incident sunlight (e.g., Vokrouhlický & Bottke 2001), but we refrain from enter this level of complexity, assuming instead that A is a constant representative value of the the sunlight reflection by the surface. We note that we also neglect the cooling effects of the sublimation processes at the surface.

The main obstacle to an analytic solution of the problem is the fourth-order, non-linear term in Eq. (3). The traditional way-around this issue is to assume that the temperature variations throughout the body are small compared to its mean value ( ). In this case, the difficult boundary term can be linearized in terms of a small parameter (quadratic and higher-order terms omitted). This approach has been carried out by a number of authors using different mathematical tools. Our description directly follows from the notation in Vokrouhlický (2006,1998,1999). For that reason, we only briefly summarize the solution referring to these indicated papers for more details.

The insolation term on the right-hand side of Eq. (3) can be expressed using a coupled series of spherical harmonics in parameters and a Fourier series in time t

 (4)

where is the solar flux at the given heliocentric distance and is the colatitude of the solar direction in the body-frame system introduced above. We recall that the diurnal approximation adopted in this paper implies that is a constant, time-independent parameter. The solar colatitude-dependent functions bnk are given by (e.g., Vokrouhlický 2006)
 b00 (5) b10 (6) (7) bnk (8)

where the last row applies to , and Pn and Pnk are Legendre polynomials and the associated Legendre functions, respectively. We note that Pn(0)=0 for n odd such that only even-degree terms appear in the insolation (4) for . Thus, the dipole term with n=1 is the only imprint of the north-south asymmetry.

Similarly to the insolation term, the temperature field can also be expressed using a mixed Fourier and spherical-harmonics expansion

 (9)

We note we are mainly interested in the temperature variations because only they can result in thermal stresses. The amplitude functions tnk(r,t) can again be given in Fourier series
 tn0(r,t) = Cn0 rn, (10) tnk(r,t) = (11)

where jn(z) denote the spherical Bessel function of order nand complex argument . We note that the zonal terms (n = 0) represent a stationary and axisymmetric field with no time dependence. The radial profile is simply given by a polynomial for zero frequency terms and jn-functions for the time-dependent terms. A simple analysis shows that the latter terms show a quasi-exponential decay in depths below the surface; is thus appropriately called the penetration depth of the diurnal thermal wave. Since , is typically much smaller than and the generic temperature field of a rapidly rotating meteoroid will have a static, but not isothermal, core, and a thin, dynamic surface layer. Only when the rotation becomes decelerated (see Sect. 3.2) may the dynamic component of the temperature field penetrate more deeply into the meteoroid.

Finally, the numerical coefficients Cn0 and Cnk in Eqs. (10) and (11) have to be determined from the boundary-condition constraints (3). One directly obtains (see also Vokrouhlický 2006,1998,1999)

 Cn0 = (12) Cnk = (13)

where , is the thermal parameter, , and

 (14)

In addition, is given by , and the average temperature . With this formulation, we retain all radial and angular information about the temperature field. Most importantly, the decomposition into zonal and non-zonal terms allows us to separately study the effects of the (quasi-)static field penetrating within the whole meteoroid and, if necessary, analyze the importance of dynamic stresses in its surface layer.

### 2.3 Thermal stress field

We now turn to the formulation of the stress field including the thermal component (see Boley & Wiener 1960; Kupradze et al. 1979; or Turcotte & Schubert (2002) for general discussion of the topic). We assume a model for the homogeneous and isotropic body undergoing small deformations characterized by the components of a symmetric strain tensor

 (15)

where is the displacement vector and index t represents the transverse operation of the appropriate matrix compound of partial derivatives of the displacement vector. While leaving our formulation in a general vectorial form in this section, we later use the system of orthonormal spherical coordinates to express components of . Using the model of a linear and isotropic solid, the relation between components of the stress tensor and the strain tensor is given by Hook's law. Thermal gradients contribute merely by a volume expansion terms, that, in the linearized theory are proportional to the temperature variation T-T0 with respect some reference value T0. Placing all the terms together, we have a generalized Hook's law in the form

 (16)

where and are the Lamé parameters from traditional linear elasticity theory and  is the linear coefficient of thermal expansion. All of these parameters have to be assumed to be dependent on the chosen reference temperature T0, which we conveniently identify with the average temperature of the body introduced in the previous section. Thus, the explicit thermal component in the stress tensor may be directly expressed using Eq. (9). There is, however, also an implicit component of the thermal contribution in the total stress field given by Eq. (16) by its contribution in the strain (deformation) field. This has to be determined by solving the Duhamel-Neumann relation (e.g., Turcotte & Schubert 2002; Kupradze et al. 1979; Boley & Wiener 1960)

 (17)

which for vanishing volumic forces  (such as the fields of solar tides or inertial forces in a rotating body) reduces to a simpler form

 (18)

In the presence of temperature gradients , is no longer a valid solution and needs to be determined by solving Eq. (18). Its uniqueness follows from the regularity of the solution in navrhuje of the whole volume of the solid and a boundary conditions (free, unload surface)

 (19)

where is again the outer normal to the surface (as in Eq. (3)).

We provide some details of the general solution in the Appendix and outline here its major steps:

• we first determine a sufficiently general form of the solution of the homogeneous form of Eq. (18);
• next we determine a particular solution of the inhomogeneous form of Eq. (18);
• because the Duhamel-Neumann Eq. (18) is linear, the general solution is a simple superposition of the two previously determined modes: ;
• we finally seek constants (A,B) such that the resulting stress field , once is substituted in Eq. (16), satisfies the boundary condition in Eq. (19).
Once we have determined the constants (A,B), we obviously also obtain a solution of the stress field throughout the whole volume of the body and can analyze whether conditions for material failure (Sect. 2.4) are reached at some point.

Because of the linearity of the Duhamel-Neumann equation and the convenience with which we can develop the thermal gradients on its right-hand side, the mixed spherical harmonics, and Fourier series, in a similar way as the temperature field itself in Sect. 2.2, we express the displacement vector in the same type of series. The only difficulty now is that we have to use vectorial spherical harmonics for instead of scalar spherical harmonics for T. This is, however, a well-known analysis and we shall use a decomposition into spheroidal and toroidal modes as discussed, e.g., in Kaula (1968) or Bullen (1975). The temperature gradient in the right-hand side of Eq. (18) is indeed a pure spheroidal field and thus we do not need the toroidal component in . The displacement vector will consist of two spheroidal terms because the differential operator on the left-hand side of the Duhamel-Neumann equation in general produces their mixture.

The Fourier (temporal) part of the series for is even simpler. Very conveniently, all zonal (axisymmetric) modes are stationary, while only the non-zonal modes represent time-dependent (periodic) part. As we saw in Sect. 2.2, the former part always penetrates throughout the whole volume of the body and thus defines the stationary stress field (only to be modified by seasonal effects), while the latter part is typically confined to a surface layer.

### 2.4 Material failure criterion

The proper goal of our paper is to seek conditions in which the thermal gradients, and the associated thermal stress field, exceed the material strength and cause a failure such as crack formation or the entire disruption of the body. Defining this critical condition accurately is a difficult task and we shall confine to a simple version known as the Griffith criterion for the brittle fracture (see, e.g., Paterson & Wong 2005). The convention used here (only for the failure criterion) is that the tension has a negative sign, whereas a positive sign represents a compression. We define and to be the maximum and the minimum eigenvalues of the stress tensor - (i.e., principal stresses) respectively. The Griffith criterion for brittle failure is then expressed by a single parametric condition (e.g., Paterson & Wong 2005)

 (20) (21)

where is the uniaxial tensile strength determined, or estimated, from experimental data. An example of the critical curve given by Eqs. (20) and (21) in the space is shown in Fig. 1. In this case we have  MPa, a value that corresponds to the estimated strength of typical ordinary chondrite material (Sect. 2.5.1).

 Figure 1: Griffith's criterion for the brittle fracture of a material with a tensile strength  MPa: and are the maximum and the minimum eigenvalues of stress tensor - . Note the positive values of  describe a compression, while the negative values describe a tension. The material (ordinary chondrite in this case) fails above, and to the left of the critical curve shown in the figure. Open with DEXTER

The three principal components of the stress tensor are computed using the Jacobi method (Press et al. 1992), which provides an efficient for determining their value at any location. We obviously use analytical expression for the stress tensor given in Sect. 3.1.

### 2.5 Material parameters

The quantitative evaluation of the thermal stresses based on the theory outlined above, and results reached below, require a number of material parameters to be known or at least estimated. Classifying them according to their relevance to the temperature or stress-field solutions, we have two groups: (i) the bulk density , the thermal conductivity K, the heat capacity c, the albedo value A, and the infrared emissivity primarily necessary for the temperature solution (Sect. 2.2); and (ii) the Lamé's parameters and , the linear coefficient of thermal expansion , and the uniaxial tensile strength  primarily necessary to determine the volumic thermal stresses (Sect. 2.3). Moreover, since several of these parameters exhibit a strong temperature dependence, we would ideally need to know them as a function of T in the range from a few hundred K to more than a thousand K relevant for the 0.05-0.2 AU heliocentric distances.

Because of the large diversity between the meteorite and meteoroid physical properties, we consider two grossly different cases (e.g., Ceplecha et al. 1998): (i) ordinary chondrite-like (OC for short) material analogous to the group I meteoroids; and (ii) primitive carbonaceous chondrite-like (CC for short) material analogous to the group III meteoroids. Some physical parameters of OC and CC materials have been determined from direct laboratory measurements of the corresponding meteorite samples (e.g., Yomogida & Matsui 1983; Britt & Consolmagno 2003; Medvedev et al. 1985). However, since measurements of different physical parameters require different instrumentation and often use different samples, no series of experiments has provided us with the complete information we need to achieve what. We therefore inspected the available literature and tried to adopt typical parameter values, and their temperature dependence, for the given class of objects, generally compiling data from several sources. These are to be considered sort of median values with a scatter of a factor of few within typically an order of magnitude.

#### 2.5.1 Ordinary chondrite (group I) material

The room-temperature measurements of , K, , , and of a H5-class meteorite Pultusk were taken from Medvedev et al. (1985). To extrapolate them to higher temperatures, we used data provided by Anderson et al. (1991), who measured them for a forstetite (taken here as an analog of the chondrite-like material). We ensured that there is a good match between the room-temperature values in the two series of measurements. The resulting relationships used in our work are

 = (22) = (23) = (24)

We then used measurements of the linear thermal expansion parameter given by Anderson et al. (1991) and preformed a linear fit approximation

 (25)

The same reference provided us with temperature-dependence measurements of the forsterite's specific heat capacity c (J K-1kg-1)

 (26)

for T<900 K, and

 c = 1020.4+0.21 T (27)

for T>900 K. The above given functional dependence c(T) yields a value 847 J K-1kg-1  for T=300 K, while Medvedev et al. (1985) obtained c=830 J K-1kg-1  for the Putulsk meteorite, indicating a good match. We assumed the thermal conductivity depends on temperature as 1/T in the studied temperature range because of interference of phonons (e.g., Opeil et al. 2010). This dependence can be described by the formula:

 (28)

Here we assumed  W m-1 K-1  is the thermal conductivity of the Putulsk meteorite for a temperature  K.

The tensile strength value, for which we did not find a direct temperature-dependence measurements, was approximated with the temperature dependence for the shear strength (Ohnaka 1995; Rocchi et al. 2004)

 (29)

where Q is the activation energy, is the gas constant, and T1 is the reference temperature. We adopt  K and T1=2500 K given by Ohnaka (1995). The constant was assumed to correspond to the measured value of the Putulsk meteorite, namely  MPa.

Finally, the characteristic albedo value A was assumed to be 0.15, typical of S-type asteroids, and the thermal emissivity was estimated to be 0.85.

#### 2.5.2 Carbonaceous chondrite (group III) material

Fewer experimental data are available for the fragile CC material. We adopted the mean bulk density for CI meteorites determined by Britt & Consolmagno (2003) and the Lame's parameters for Axtell CV3 meteorite determined by Flynn (2004). Owing to the lack of data, we assumed the same temperature dependence of these quantities as given above for the forsterite (Anderson et al. 1991)

 = (30) = (31) = (32)

We were unable to find an appropriate laboratory measurement for the other parameters we need for this class of objects. In this situation, we basically adopted the same parametric relationships as given above for the OC class with the following modifications: (i) the nominal tensile strength in Eq. (29) was given a value of -2 MPa to express the weaker nature of CCs compared with OCs (the assumed order of magnitude difference between the values roughly corresponds to the difference between the observationally-determined dynamic pressures of the group I and III meteoroids; e.g., Ceplecha et al. 1998; Borovicka 2007); (ii) the thermal conductivity  W m-1 K-1 of the Cold Bokkeveld CM2 meteorite for temperature  K (Opeil et al. 2010) was considered in Eq. (28); and (iii) the albedo value A=0 and thermal emissivity was considered in this case to be consistent with the overall darker character of the CC material.

## 3 Results

### 3.1 Expressions for the thermal stress tensor

We summarize our final expressions for the thermal stress tensor components based on the formulation given above and the intermediate results given in the Appendix. Since we use the system of spherical coordinates to parameterize the dependence of both temperature T and the stress tensor on position, it appears most efficient to project the components of the latter onto the system of three orthonormal vectors

 (33)

where is the position vector of a chosen point in the body. We thus compute the components , and so on. Since the outward normal to a sphere is equal to , the boundary conditions in Eq. (19) are easily expressed by the three relations .

As noted above, the thermal deformation, expressed by the displacement vector , is given in terms of a mixed vectorial spherical harmonic series in coordinates and a Fourier series to represent its time dependence. Applying the differential operators in Eqs. (15) and (16) that relate the displacement vector to the stress tensor, we obtain components of the resulting stress tensor in terms of a series of spherical functions and their derivatives. Conveniently, the structure of the Fourier series in time remains the same. In particular, all axisymmetric (zonal) terms are stationary, whereas all non-axisymmetric terms are time-dependent and correspond to a frequency for the order k in the spherical harmonics functions.

#### 3.1.1 Time-independent part of the thermal stress tensor

We begin with the stationary part of the solution since it allows more compact final expressions for the components of the stress tensor and this is often a more important part of the deformation.

We introduce auxiliary functions

 (34)

where x=r/R and

 (35)

we suppressed the parametric dependence of the radial profile G-functions keeping note only of its r-dependence. With these functions defined, we then have
 = (36) = (37) = (38) = (39)

All other components of the stress tensor vanish. We note that the -dependence in the argument of the spherical functions above is dummy because here we always have zonal terms only and thus the field functions are truly axisymmetric. A few comments are in order.

We note the above given series start with the quadrupole term (n = 2) only, while for instance the temperature field in Eq. (9) contains also the dipole part. This is because the corresponding linear temperature field across the body produces a deformation equivalent to a linear displacement that is not capable of causeing stress (e.g., Iesan & Scalia 1996). The bulk of the body is mainly affected by the quadrupole part because higher multipoles have gradually steeper xn-2 decays for small r. We should also point out that because its amplitude G2 is proportional to , the quadrupole contribution vanishes near the node of the second-degree Legendre polynomial. This effect produces anomalously small stresses in the body when the solar direction is tilted by this angle with respect to the spin axis of the body. The boundary conditions are readily satisfied by the (1-x2)term in and ( is nil). In terms of physical parameters, the stress field is (i) linearly proportional to the thermal expansion parameter ; (ii) roughly linearly proportional to the Lame's parameters and (both of which have the same order of magnitude); and (iii) increases with increasing R in some range of values, near , up to a saturation value. The last property reflects that , and thus it depends on R only through the R-dependence of in the denominator factor in Eq. (12). Because , the thermal stresses in Eqs. (36) to (39) increase roughly inversely proportionally to the square root of the heliocentric distance making them more important close to the Sun. This is obviously a fairly expected trend.

The stationary part of the stress field readily dominates in the two limits: (i) very fast rotation of the body ( ); and (ii) solar direction along the spin axis of the body ( or ). In the first case, the penetration depth of the diurnal thermal wave shrinks to zero and all time-dependent components of the field are pushed to an infinitesimal slab near the surface of the body. In the second case, the axisymmetry is clearly imposed by the geometry of the problem. In terms of a mathematical description, we note that the time-dependent part of the field is proportional to the coefficients that vanish for and .

#### 3.1.2 Time-dependent part of the thermal stress tensor

We now consider the time-dependent part of the deformation and the stress field, which is expressed by the non-zonal () terms. In this case, we were unable to derive as compact a form of the final result as for the stationary part and we give it in a piecewise form.

We mentioned above, in Sect. 2.3, that the total stress tensor is composed of two parts

 (40)

the first one of which, , depends on the deformation field in the usual way given by the Hook's law, and the second part being the explicit temperature-dependent term in the generalized Hook's law (16). The deformation vector is a linear superposition of two free, spheroidal modes and , derived in Sect. A.1, and a particular solution of the complete Duhamel-Neumann equation, given in Sect. A.2. All of them share the spherical-harmonics development and Fourier-development structure introduced above. We may therefore write

 (41)

for the harmonics term of degree n, order k, and frequency k in the Fourier space. The constants Q1nk and Q2nk need to be determined to satisfy the surface boundary conditions (19) and Sect. A.3. We then spectrally decompose the stress-tensor field

 (42)

where

 (43)

The explicitly temperature-dependent stress field is given by , and its spectral decomposition follows directly from the decomposition of the temperature field given in Sect. 2.2.

In the next few paragraphs, we provide expressions for the individual terms in Eq. (42). We list the non-zero components of the stress tensor field only and provide their projections onto the orthonormal basis given by Eq. (33). To make the notation shorter, we drop the index nk of the stress tensor field components, which should nevertheless be understood from their presence in right-hand sides of the formulae given below.

First, we consider the components of and corresponding to the two spheroidal modes that provide a solution to the homogeneous Duhamel-Neumann equation. For the first mode, we obtained

 = (44) = (45) = (46) = (47) = (48) = (49)

 Figure 2:Thermal stress and temperature distribution in a meridional section of a centimeter-size CC meteoroid at 0.14 AU heliocentric distance. The Sun (indicated by the arrow) is located above the north pole; as a result, the time-dependent component of both fields vanishes. Isolines show given constant values of the ( left panel) and ( middle panel) principal components of the stress tensor, both in MPa, and the temperature in K ( right panel). Open with DEXTER
The computations are simpler for the divergence-free second spheroidal mode ( ), for which we obtain
 = (50) = (51) = (52) = (53) = (54) = (55)

We next indicate the thermal stress components of the particular solution of the Duhamel-Neumann equation: . We find it suitable to define

 (56)

 (57)

With them, we have
 = (58) = (59) = (60) = (61) = (62) = (63)

Finally, we recall that the explicit temperature contribution in the generalized Hook's law is a pure volumic expansion

 (64)

Equipped with these partial results, we may now compute the coefficients Q1nk and Q2nk. The task is separate for the dipole terms (n=1) and the remaining multipoles (). This is because the surface components of the second spheroidal mode vanish, i.e., , for the dipole term.

In the first case, we find that the surface condition expressed in terms of the rr, r, and r components of the stress tensor are all linearly dependent. For this reason, they represent a single condition that yields

 (65)

The coefficients are dummy and disappear from the final expression for the dipole part of the thermal stresses because all components of are zero: (i) the rr, r, and r are readily proportional to (n-1)=0; while (ii) the , , and components vanish because of the identities of the spherical functions.

For the higher multipole terms, the rr and r surface conditions are linearly independent and provide a constraint on the Q1nk and Q2nk constants. After solving the corresponding set of algebraic equations, we obtain

 Q1nk = (66) Q2nk = (67)

### 3.2 Examples of thermal stresses in small meteoroids

Before investigating the parameter dependencies of the conditions for the onset of thermal fission of small meteoroids, we illustrate the thermal-stress fields in meteoroids in several individual cases. Obviously, the previous analytical formulation is used and evaluated using a numerical implementation in Fortran. We choose a CC meteoroid of 1 cm diameter at a distance of 0.14 AU from the Sun. Material constants are those from Sects. 2.5.1 and 2.5.2.

 Figure 3:Thermal stress and temperature distribution in a meridional section of a centimeter-size CC meteoroid at 0.14 AU heliocentric distance. The Sun (indicated by the arrow) is located above the equator) and the rotation frequency is assumed to be 100 Hz. Isolines show given constant values of the (left panel) and (middle panel) principal components of the stress tensor, both in MPa, and the temperature in K (right panel). Open with DEXTER

We first assume a polar direction to the Sun, i.e., . In this case, the time-dependent part of the stress field vanishes and we are left with only the stationary and axisymmetric part. The distribution of the maximum and minimum principal values of the stress tensor (see Sect. 2.4) in a meridional section is shown in Fig. 2. The extreme values of and are found on the surface and at the equator (8 MPa and -20 MPa, respectively). Moreover, has a local extreme at the center of the body (-7 MPa). The temperature distribution (right panel of Fig. 2) ranges from 1140 K at the north pole to 580 K at the south pole.

We next consider the situation in which the Sun is above the equator ( ) and the meteoroid's rotation frequency is 100 Hz. In this case, both stationary and time-dependent stress fields exist, but the latter is confined to a thin surface layer of an approximate width  mm for this high rotation frequency (recall that ). The stationary component of the stress field thus again dominates across the major part of the meteoroid's volume. Figure 3 shows the distribution of the maximum and minimum  principal values of the stress tensor in a meridional section containing the solar-direction vector. The eigenvalue has a maximum of 6 MPa near the surface at the equator and a local maximum at the center of the body (4 MPa). The eigenvalue  has extreme values at poles (-6 MPa) and a local maximum in the equatorial plane. The assumed very high rotation frequency makes the temperature field nearly axisymmetric with minimum values 660 K at the poles and maximum values 770 K at the equator.

 Figure 4: Thermal stress and temperature distribution for a slowly rotating centimeter-size CC meteoroid at 0.14 AU heliocentric distance. The Sun (the direction of which is indicated by the arrow in the left panels) is located at the colatitude with respect the spin axis (oriented upwards) and the rotation frequency is assumed to be 0.1 Hz. Left panels show a meridional section of the body (containing the solar direction), right panels show an Aitoff projection of the surface. Isolines show given constant values of the (top panels) and (middle panels) principal components of the stress tensor, both in MPa, and the temperature in K (bottom panels). Open with DEXTER

Perhaps the most interesting and complex case arises for slowly rotating meteoroids. In our last example, we thus assume the rotation frequency of 0.1 Hz only and that the Sun is at a colatitude, both for a 1 cm CC meteoroid at a 0.14 AU heliocentric distance. The chosen solar colatitude minimizes the role of the stationary component of the stress field (see Sect. 3.1.1), while the slow rotation now allows the time-dependent component of the stress field to penetrate more deeply into the body. Figure 4 shows that the surface of the insolated hemisphere is loaded by compression (the maximum value of  MPa is close to the subsolar point) and the surface of the shadowed hemisphere is loaded by tension (extreme value of is almost -30 MPa). The major variations in the temperature and the stress field penetrate approximately to the depth of 1 mm. The maximum temperature 900 K occurs near the subsolar point, only a minor longitudinal shift being caused by the finite value of the thermal conductivity (bottom and right panel of Fig. 4), while the minimum temperature is on the south pole (600 K). The amplitude of temperature variations at the subsolar latitude during one rotation cycle is 140 K.

### 3.3 Destruction of small meteoroids

Except for the material properties of the meteoroid, which we keep fixed and the same as in Sect. 2.5, there are four major parameters in our model that determine the thermal stress magnitude and therefore also the possibility of thermal fission: the heliocentric distance a (translating into the reference temparature value), the size of the body D (or its radius R), the rotation frequency f (or ), and the instant solar colatitude . We assume that cracks and fissures start to form and propagate from a point in the body where the Griffith failure criterion (Sect. 2.4) is satisfied. This process very likely has its onset near the surface of the body and may have a complex influence on how the thermal fission process propagates further into the body. This is because the fractured surface layer may start to thermally shield the interior of the body by being of lower thermal conductivity. This dynamical model of fission is, however, beyond the scope of this paper and we plan to face some of its aspects in a forthcoming publication. In this work, we take a simpler standpoint and assume that if the material failure conditions take place at the center of the body, the conditions for fracturing are also fulfilled in the majority of its volume and a catastrophic thermal burst occurs. In the following sections, we thus compute the principal stress-tensor components and at the center of meteoroid for the two extreme composition cases of OCs and CCs. Most often, it is the value (thermal tension) that overrides the Griffith-criterion line (Fig. 1).

Nevertheless, we are left with the four-dimensional dependence on a, D, f, and parameters. To illustrate the effects in a simple way, we fixed two of these parameters and changed the other two in the next few sections.

 Figure 5: The principal component - of the stress tensor at the center of a meteoroid as a function its size and the heliocentric distance. The isobars are labeled according to their values in MPa. The upper plots for the solar colatitude and bottom ones for . The left panels for the CC material parameters, the right panels for the OC material parameters. The rotation frequency f is inversely proportional to the size D as in Eq. (68). The shaded zone indicates where the critical line of the Griffith criterion has been overrun and the thermal fission occurs. Open with DEXTER

#### 3.3.1 Dependence on the heliocentric distance and size

We first analyze the onset of the thermal fission in the plane of heliocentric distance versus size parameters. We choose two values of the solar colatitude, and . While in the first case the results do not depend on the rotation frequency, in the second case we have to assume some specific value of f. Instead of using one particular value of f, we consider instead a parametric dependence f(D)such as

 (68)

where D is in meters and f in Hertz. This relation of inverse-proportionality matches very roughly suggested rotation periods of centimeter- to decimeter-size meteoroids reported in
Beech & Brown (2000) and extends to Ceplecha's (1996) determination of the 3.3 s rotation period of the Lost City fireball. We consider sizes of between 1 mm to 10 m, again recalling that our formulation, which disregards seasonal thermal effects, prevents us from considering objects larger than a few metres in size. Heliocentric distances range from 0.05 AU to 1 AU in our example. Figure 5 shows the values that we obtained at the center of the CC- and OC-type meteoroids. Overall, increases with the increasing size of the body and with decreasing heliocentric distance (note that the material strength also decreases with decreasing heliocentric distance because of its temperature dependence). The principal take-away message here is that, especially for weak material such as CCs, a size limit exists at a given heliocentric distance (0.3 AU, say) above which the meteoroids get promptly fissioned (see further discussion in Sect. 4).

A more detailed comparison indicates that the case (upper panels in Fig. 5) caused more significant thermal stresses than to the case (bottom panels in Fig. 5). This is because at non-zero values the surface layer absorbs some of the thermal gradients and prevents the building of a large global gradient across the whole body.

 Figure 6: The principal component - of the stress tensor at the center of a meteoroid as a function its size and the solar colatitude . The isobars are labeled by values in MPa. The left panel for the CC material parameters, the right panel for the OC material parameters. The rotation frequency f is inversely proportional to the size D as in Eq. (68) and the heliocentric distance is 0.14 AU. The shaded zone indicates where the critical line of the Griffith criterion has been overrun and the thermal fission occurs. Open with DEXTER

#### 3.3.2 Dependence on the solar colatitude and size

We next consider the dependence of the thermal-stress field on the solar colatitude , and the size D of the meteoroid, fixing the heliocentric distance to be 0.14 AU (perihelion distance of the Geminid stream) and assuming the f(D) relation from Eq. (68). Sizes are in the millimeter to meter range as above and allows us to span the interval of to values. Because we consider thermal stresses at the center of the body that are generally unreached by the time-dependent component of the stress field, which is confined to the surface layer for our f(D)values, there is a symmetry between and results.

Figure 6 shows the results. From the analytic analysis in Sect. 3.1.1, we can infer that the stress field basically vanishes for . More importantly, there is a general trend toward more relaxed stress field when going from the configuration with the Sun above the rotation pole ( ) to the configuration with the Sun above the equator. As mentioned above, larger values imply that some of the thermal gradients are absorbed by the surface (dynamic) layer with less damage occurring in the bulk of the body.

 Figure 7: The principal component of the stress tensor in the center of a meteoroid as a function its size and the rotation frequency f. The isobars are labeled by values in MPa. The left panel for the CC material parameters, the right panel for the OC material parameters. The heliocentric distance is 0.14 AU and the chosen value of the solar colatitude maximizes the thermal stresses (Fig. 6). The solid line shows the f(D) relation from Eq. (68). The dashed line indicates the situation in which the radius of the body is five times the penetration depth of the diurnal thermal wave (Eq. (69)). The shaded zone indicates where the critical line of the Griffith criterion has been overrun and the thermal fission occurs. Open with DEXTER

#### 3.3.3 Dependence on the rotation frequency and size

Finally, the dependence of the thermal stresses at the center of the meteoroid on the rotation frequency f and the size D can be seen in Fig. 7 (we now relax the strict f(D) relation from Eq. (68), which is shown by the solid line). For illustration, we assumed the solar colatitude and the heliocentric distance a=0.14 AU.

We note a general pattern for a range of different sizes and the two material classes (CCs and OCs): is almost constant with decreasing f until a transition to a higher, but also constant value, occurs for f lower than some critical value. One can readily explain this behavior by recalling that at high frequencies, the time-dependent part of the stress field is confined to a tiny surface layer and the frequency-independent stationary part of the field determines the stress field at the center. When the frequency decreases below a critical value for which , thus when the diurnal thermal wave reaches the center of the body, it is basically the time-dependent component of the stress field that dominates the effect. Its peak value during one cycle is roughly independent of the already low rotation frequency. The transition between the two regimes is thus expressed by

 (69)

shown by the dashed line in Fig. 7. This situation corresponds to when the radius of the meteoroid is equal to five times the penetration depth of the diurnal thermal wave, indicating that the time-dependent stress field reaches somewhat below . We note the empirical (observed) f(D) relation from Eq. (68) is moved toward the range of smaller stresses for the given size. In other words, the observed meteoroids tend to rotate more rapidly and prevent potentially large thermal stresses in their volume by not allowing the diurnal thermal wave to penetrate sufficiently inside the body.

## 4 Conclusions and further work

We have used a simple analytical model to evaluate thermal stresses in small meteoroids near the Sun. The most notable approximations are (i) a spherical shape with a homogeneous distribution of physical parameters in the body; and (ii) regular rotation about the spin axis fixed in space. We also did not include seasonal thermal effects in our model, thus the meteoroid's revolution about the Sun. As a result, only bodies smaller than few metres are eligible for our theory. It should be, however, noted that the absence of the seasonal effects is the least difficult extension of our theory. To zeroth order, it would suffice to modify the zonal (k=0) insolation terms in Eq. (4) to ensure time-variability with the mean motion frequency, and express them along the lines given in Eqs. (8) and (9) of Vokrouhlický (1999). Conceptually more interesting efforts should perhaps be directed towards removing the (i) and (ii) simplifications above. We plan to do so in a forthcoming paper.

Given the lack of these modeled generalizing features, such as the potential effects of insulating (granular) layer on the surface of the meteoroid, we also hesitate to directly compare our results with the observed facts. We nevertheless provide some initial discussion below, while more quantitative analyses are postponed to a future paper.

The major results obtained above are twofold. First, meteoroids at small heliocentric distances may develop large enough thermal stresses to induce fission above some size limit. For instance, our results for a conservative case of CC material with shown in Fig. 5c would imply that 10 cm size meteorids would thermally fission below 0.25 AU heliocentric distances. Second, observed meteoroids tend to occupy a safe suite in the rotation frequency versus size space by maintaining relatively rapid rotation. However, if some processes, such radiation torques known as the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect (e.g., Capek & Vokrouhlický 2004), were able to decelerate their rotation they would become more vulnerable to the thermal fission. This is because the diurnal thermal wave would penetrate more deeply into the body and carry along steeper thermal gradients affecting thus most of its volume. In the same time, rapid rotation of the meteoroids keeps the diurnal thermal wave with its associated large thermal gradients very near the surface possibly resulting in a fractured layer of lower thermal conductivity. This insulating shell may in term protect the bulk of the meteoroid from large thermal stresses and thus play a self-regulating role. The question then is from which size the meteoroids are able to build and retain this protective layer. We also recall that while we have assumed some reasonable values for the thermal and strength parameters, an uncertainity of a factor of a few is quite possible. For instance, if the CC's critical tensile strength was twice as large as our assumed value, meteoroids of all sizes would survive approaches to 0.05 AU in the example shown in Fig. 5c. What do we know about meteoroids in the decimeter- to dekameter- size range from the observations and how do they compare with our theoretical results? (We focus in this respect, on the weak objects with cometary origins.)

Several CC meteorite falls have been well documented and interpreted. Exceptional example of these events are the cases of Tagish Lake (e.g., Brown et al. 2000) and Orgueil (e.g., Gounelle et al. 2006) meteorites. In both cases, the estimated parent object sizes are close to dekameter values and both pre-atmospheric orbits have perihelia of 0.87 AU. Rotation periods of several seconds, as in the Lost City meteorite case (Ceplecha 1996), and not in the vicinity of would protect these bodies against thermal disruption even without assuming the insulating (granular) surface layer. We would need to assume it only for perihelia below 0.3-0.4 AU in these cases.

In addition, a large amount of data about sporadic meteors is available from the automated fireball networks. As an example, we mention the results from the European networks overviewed by Oberst et al. (1998). Focusing on the brightest recorded bolides, we note several group-IIIb meteors of an estimated meter-range pre-atmospheric size and low perihelia. For instance, Visla bolide had a perihelion at 0.22 AU. Inspecting results in Fig. 5, we conclude that this object should have had a larger tensile strength than we assumed for the CCs meteorite class and/or was protected from thermal bursting at the perihelion of its pre-atmospheric orbit by an existing insulating layer on the surface. Certainly more analysis is needed to exploit the wealth of the network data of sporadic meteors, but we postpone this to the future paper after complementing the current work with the effects of aninsulating surface layer and/or the seasonal-thermal wave damping into the volume of the meteoroid.

In the future, a large amount of information could potentially come from observations of meteoroids in streams and documented cometary fragmentations. The regular observation of stream meteoroids, as they enter and disintegrate in the atmosphere, typically provides evidence of particles smaller than between centimeters and decimeters. Searches for larger components in these streams were reviewed by Beech & Nikolova (2001). For instance, Perseid and Leonid streams, both associated with long-period comets, may show observational evidence of meter- to dekameter- sized fragments. With their perihelia at 1 AU, we can still ensure that they withstand thermal stresses, especially if the surface is covered with a thin, insulating layer.

Some other meteoroid streams have smaller perihelion distances, such as -Aquarids 0.07 AU, Geminids 0.14 AU, and/or Monocerotids 0.19 AU. We are not aware of information about the largest-observed meteoroids in these streams, but they would be the prime targets in investigating the possible depletion by larger meteoroids in the millimeter- to decimeter- size range (especially the fragile and slowly rotating ones).

Evaluation of the subsurface thermal stress field for fragments of disrupted comets in the Kreutz family with very low perihelia might be another interesting extension of our work. For that project, we would need to include modeling of the seasonal effects along a very eccentric orbit. Semi-numerical methods may be necessary to determine the amplitude of the different, and numerous, Fourier terms of the incident solar flux, but once these are established, results from this paper may be used to evaluate the thermal stress field.

Acknowledgements
We thank Guy Consolmagno for comments on the first version of this paper. The work of D.C. was supported by the Grant Agency of the Czech Republic under a contract 205/09/P455. The work of D.V. was supported by the Research Program MSM0021620860 of the Czech Ministry of Education.

## References

1. Anderson, O. L., Isaak, D. L., & Oda, H. 1991, J. Geophys. Res., 96, 18037 [NASA ADS] [CrossRef] [Google Scholar]
2. Beech, M., & Brown, P. 2000, Planet. Space Sci., 48, 925 [NASA ADS] [CrossRef] [Google Scholar]
3. Beech, M., & Nikolova, S. 2001, Planet. Space Sci., 49, 23 [NASA ADS] [CrossRef] [Google Scholar]
4. Boley, B. A., & Wiener, J. H. 1960, Theory of thermal stresses (New York: Wiley) [Google Scholar]
5. Borovicka, J. 2007, in IAU Symp., ed. A. Milani, G. B. Valsecchi, & D. Vokrouhlický, 236, 107 [Google Scholar]
6. Britt, D. T. & Consolmagno, G. J. 2003, Meteoritics & Planetary Science, 38, 1161 [CrossRef] [Google Scholar]
7. Brown, P. G., Hildebrand, A. R., Zolensky, M. E., et al. 2000, Science, 290, 320 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
8. Bullen, K. E. 1975, The Earth's density (London: Chapman and Hall) [Google Scholar]
10. Ceplecha, Z., Borovicka, J., Elford, W. G., et al. 1998, Space Sci. Rev., 84, 327 [NASA ADS] [CrossRef] [Google Scholar]
11. Capek, D., & Vokrouhlický, D. 2004, Icarus, 172, 526 [NASA ADS] [CrossRef] [Google Scholar]
12. Capek, D., & Borovicka, J. 2009, Icarus, 202, 361 [NASA ADS] [CrossRef] [Google Scholar]
13. Flynn, G. J. 2004, Earth, Moon, and Planets, 95, 361 [Google Scholar]
14. Gounelle, M., Spurný, P., & Bland, P. A. 2006, Meteoritics & Planetary Science, 41, 135 [CrossRef] [Google Scholar]
15. Grinin, V., Natta, A., & Tambovtseva, L. 1996, A&A, 313, 857 [NASA ADS] [Google Scholar]
16. Hall, K. 1999, Geomorphology, 31, 47 [NASA ADS] [CrossRef] [Google Scholar]
17. Iesan, D., & Scalia, A. 1996, Thermoelastic Deformations (Dordrecht: Kluwer Academic Publishers) [Google Scholar]
18. Jenniskens, P. 2006, Meteor Showers and their Parent Comets (Cambridge University Press, Cambridge) [Google Scholar]
19. Kaula, W. M. 1968, An introduction to planetary physics - The terrestrial planets (New York: Wiley) [Google Scholar]
20. Kuehrt, E. 1984, Icarus, 60, 512 [NASA ADS] [CrossRef] [Google Scholar]
21. Kupradze, V. D., Gegelia, T. G., Basheleishvili, M. O., & Burchuladze, T. V. 1979, Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity (New York: North-Holland Pub. Co.) [Google Scholar]
22. Marchi, S., Delbó, M., Morbidelli, A., Paolicchi, P., & Lazzarin, M. 2009, MNRAS, 400, 147 [NASA ADS] [CrossRef] [Google Scholar]
23. Medvedev, R. V., Gorbatsevich, F. I., & Zotkin, I. T. 1985, Meteoritika, 44, 105 [NASA ADS] [Google Scholar]
24. Oberst, J., Molau, S., Heinlein, D., et al. 1998, Meteoritics & Planetary Science, 33, 49 [Google Scholar]
25. Ohnaka, M. 1995, Geophys. Res. Lett., 22, 25 [NASA ADS] [CrossRef] [Google Scholar]
26. Opeil, C. P., Consolmagno, G. J., & Britt, D. T. 2010, Icarus, 208, 449 [NASA ADS] [CrossRef] [Google Scholar]
27. Parkus, H. 1976, Thermoelasticity (Wien - New York: Springer-Verlag) [Google Scholar]
28. Paterson, M. S., & Wong, T. 2005, Experimental Rock Deformation - The Brittle Field (Berlin: Springer) [Google Scholar]
29. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes, 2nd edn (New York: Cambridge University Press) [Google Scholar]
30. Rocchi, V., Sammonds, P. R., & Kilburn, C. R. J. 2004, Journal of Volcanology and Geothermal Research, 132, 137 [NASA ADS] [CrossRef] [Google Scholar]
31. Shestakova, L. I., & Tambovtseva, L. V. 1997, Earth, Moon, and Planets, 76, 19 [Google Scholar]
32. Tambovtseva, L. V., & Shestakova, L. I. 1999, Planet. Space Sci., 47, 319 [NASA ADS] [CrossRef] [Google Scholar]
33. Tauber, F., & Kuhrt, E. 1987, Icarus, 69, 83 [NASA ADS] [CrossRef] [Google Scholar]
34. Thorne, K. S. 1980, Rev. Mod. Phys., 52, 299 [NASA ADS] [CrossRef] [Google Scholar]
35. Turcotte, D. L., & Schubert, G. 2002, Geodynamics (Cambridge: Cambridge University Press) [Google Scholar]
36. Urquhart, M. L., & Jakosky, B. M. 1997, J. Geophys. Res., 102, 10959 [NASA ADS] [CrossRef] [Google Scholar]
39. Vokrouhlický, D., & Bottke, Jr., W. F. 2001, A&A, 371, 350 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
40. Vokrouhlický, D. 2006, A&A, 459, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
41. Vokrouhlický, D., Nesvorný, D., Dones, L., & Bottke, W. F. 2007, A&A, 471, 717 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
42. Yomogida, K., & Matsui, T. 1983, J. Geophys. Res., 88, 9513 [Google Scholar]

## Appendix A: Solution of the Duhamel-Neumann equation

We outline the main steps needed to solve the Duhamel-Neumann equation in Eq. (18) for our work. The temperature field, which produces the thermal stresses, is assumed to have a linearized form where is given by Eq. (9). The uniqueness of the solution arises from (i) the regularity in the whole volume; and (ii) matching the free boundary conditions (given by Eq. (19)) at the surface r = R.

There are different ways in which we can decompose the displacement vector into spherical-harmonics-type expansion (see, e.g., Thorne 1980, for an insightful review). Here we use a decomposition into spheroidal and toroidal components traditionally used in geophysical analyses (e.g., Kaula 1968; Bullen 1975). With this approach, related to what Thorne (1980) calls pure-spin vector harmonics, we have

 (A.1)

with
 = (A.2) = (A.3)

The first component in Eq. (A.1), , is actually two separate spheroidal terms, that we call S1 and S2 in Sect. 3.1.2, characterized with radial-profile amplitudes Unk(r) and Vnk(r). The second component in Eq. (A.1), , is the toroidal component.

The spheroidal character of our source (temperature) term the Duhamel-Neumann equation implies two simplifications. First, the toroidal part of the displacement vector becomes negligible and we have Wnk=0. Second, we can restrict the summation over degrees n in Eq. (A.1) to the dipole and higher-order terms only, ignoring the monopole n=0. This is because the monopole part would correspond to purely radial temperature field, such as has been considered, for instance, in the previous works on our topic (e.g., Tambovtseva & Shestakova 1999; Shestakova & Tambovtseva 1997; Kuehrt 1984). Our temperature representation does not contain a non-trivial, purely radial profile and the only viable free monopole term must have U00=W00=0to match the boundary conditions. Finally, we note that we also anticipated the Fourier-development structure in Eq. (A.1) as it follows from the source ( development).

Substituting the spheroidal-vector representation of into the Duhamel-Neumann Eq. (18) we obtain the following system of equations for the radial profile of the amplitude functions Unk(r) and Vnk(r) ():

 (A.4)

and

 (A.5)

Here we found it useful to separate the r- and t-dependences of the tnk(r,t) amplitudes of the development in Eq. (9) and introduce pure radial parts Tnk(r) such that . These represent source terms in Eqs. (A.4) and (A.5).

While the solution of Unk(r) and Vnk(r) is coupled by means of Eqs. (A.4) and (A.5), the fundamental implication of the Duhamel-Neumann equation linearity is that amplitude terms of different degrees and orders in the spherical harmonics development as well as the different Fourier modes are not mixed and can be solved separately.

Once we obtain Unk(r) and Vnk(r), we can readily compute components of the corresponding stress tensor arising from the displacement vector field by using the Hook's law in Eq. (16). Given its linearity, we thus again have

 (A.6)

where . Projecting components of the stress tensor onto the orthonormal basis from Eq. (33), as outlined in Sect. 3.1, we obtain (for simplicity we dropped here the degree- and order-indexes n and k)
 = (A.7) = (A.8) = (A.9) = (A.10) = (A.11) = (A.12)

The partial derivatives of the spherical functions Ynk are computed using
 = (A.13) = (A.14)

Equations (A.7) to (A.12) yield components of the stress tensor that explicitly depend on the displacement vector . The part , which explicitly depends on the temperature (see the generalized Hook's law in Eq. (16), should be added separately to the total stress-tensor field. Because of the explicit analytical solution for , this is achieved at no computational expense.

Because of the linearity of the Duhamel-Neumann equation, a general solution is expressed in terms of a linear superposition of (i) a solution of the homogeneous system; and (ii) a particular solution of the inhomogeneous system. The next two sections discuss the two cases separately.

### A.1 Solution of the homogeneous Duhamel-Neumann equation

Equations (A.4) and (A.5) with zero right-hand sides represent the homogeneous Duhamel-Neumann equation broken into parts corresponding to the individual spheroidal modes. Its solution is quite complicated, but may be significantly simplified in our case. This is because for the range of material parameters, sizes and rotation frequencies that apply for meteoroids we always have . With these we may neglect the troublesome term in Eqs. (A.4) and (A.5). A major implication of this is then that the system of solutions of the homogeneous Duhamel-Neumann equation become degenerate in the k (order) index of the spherical-harmonics representation.

Adopting the aforementioned approximation, the homogeneous system of Eqs. (A.4) and (A.5) now has a form of Euler equations. As such, it has a fundamental system of power-law solutions Uink=Qi rmi and Vink = rmi with , real-valued exponents mi and amplitudes Qi. After a straightforward algebra, we obtain

 U1nk(r) = (A.15) V1nk(r) = rn+1, (A.16) U2nk(r) = (A.17) U3nk(r) = (A.18) U4nk(r) = (A.19) V4nk(r) = r-n. (A.20)

The last two modes, 3 and 4, diverge at the center r=0 and therefore must be excluded. We are thus left with the first two modes, 1 and 2, that produce the spheroidal modes and in Sect. 3.1.2 and whose associated stress field was given in Eqs. (44)-(49) and Eqs. (50)-(55). Obviously, they have also been used to obtain the stationary part of the stress field discussed in Sect. 3.1.1. We note that the second spheroidal mode represents a pure shear with no volumic changes (compression or expansion) because .

### A.2 Particular solution of the Duhamel-Neumann equation

We next find a particular solution of the inhomogeneous Duhamel-Neumann equation with the thermal source . We divide this task into a discussion of the stationary case (k=0) and time-dependent case (). In both cases, we again use the approximation of neglecting the terms in Eqs. (A.4) and (A.5).

#### A.2.1 Time-independent part

The stationary temperature field is given by (Eqs. (9) and (10)) and thus . We again search the fundamental system of solutions in a power-law form and with some real-valued exponents and amplitudes . After a brief algebraic derivation, we obtain

 = (A.21) = (A.22)

We note that this mode has the same radial profile as the spheroidal model found above.

#### A.2.2 Time-dependent part

The time-dependent temperature field is given by with (Eqs. (9) and (11)) and thus . We assume that the particular solution has a form . Substituting this ansätz to the Duhamel-Neumann equation Eq. (18), we obtain

 (A.23)

where we have suitably assumed that the arbitrary constant on the right-hand side canceled the monopole (constant) temperature part. This is an inhomogeneous wave equation on a sphere that, however, takes a simple form because of the spherical-harmonic and Fourier structure of the source term on the right-hand side. Assuming thus a separation

 (A.24)

we obtain

 (A.25)

and

 (A.26)

where is the elastic P-wave velocity as above. Adopting again the approximation , we may neglect the first term in the denominator of Eq. (A.26). Translating this solution into the amplitude-functions of the spheroidal-field representation (A.2), we finally obtain
 = (A.27) = (A.28)

The corresponding stress tensor is expressed by Eqs. (58)-(63).

### A.3 Complete expression of the thermal stress tensor

The complete solution of the Duhamel-Neumann equation is a linear combination of the free-spheroidal modes and from Sect. A.1 and the particular mode from Sect. A.2. In the individual spherical harmonics modes, we have , where Q1nk and Q2nk are some coefficients. We have to choose them to satisfy the surface boundary condition (19), namely at r=R. Here the total stress tensor is given by

 (A.29)

or again in the spherical harmonics modes

 (A.30)

The truly active and independent conditions are and , from which the two constants Q1nk and Q2nk follow. One easily checks that the third condition, , is always linearly dependent on (Eqs. (A.8) and (A.9)) and thus we do not need to consider it.

We were able to carry out all necessary algedraic manipulations and obtain a close form of the resulting formulae for the case of the stationary (zonal, k=0) part of the stress field. These are given in Eqs. (36)-(39) (Sect. 3.1.1). In the case of the time-dependent part of the stress field, the algebra is more involved and we could not reach as simple and compact results as for the time-independent part. We thus confine ourselves to provide formulae for the stress-tensor components of the individual components and those for the integration constants Q1nk and Q2nk(Sect. 3.1.2).

## Footnotes

... meteoroids
Appendix A is only available in electronic form at http://www.aanda.org
... temperature
We later consider their dependence on the mean temperature of the body to compare results across a range of different heliocentric distance but do not take into account their temperature dependence when solving the heat diffusion problem.
... above
Note that we always work in the reference system attached to the body such that the solar longitude is represented by the rotation-frequency-dependent term in Eq. (4); we also denote . The last two terms in Eq. (4) may also be rewritten to illustrate this property.
... dependence
It is obviously mainly in these terms that the seasonal, mean-motion variations would appear, if we were to include also the meteoroid revolution about the Sun.
... tensor
Eigenvalues of tensor - instead of are computed and this assures the sign convention mentioned in the text. Typically, is positive (i.e., compression) and negative (i.e., tension).
... match
Interestingly, the adopted c(T) dependence corresponds quite well to the model of lunar soil given by Urquhart & Jakosky (1997).
... part
This dipole part is most responsible for the overall temperature decrease from the solar-illuminated hemisphere to the opposite side of the body.
... body
The applicability of our theory obviously reduces when becomes comparable to or smaller than the granularity level of the surface.
... define
The Bessel equation provides us with an expression for the function using the previously introduced variables .
...tex2html_comment_mark
In Sects. 3.2 and 3.3, we used rotation frequency f, cycles per second or Hz, that is related to using the relation: .
... profile
Note this does not mean there is not an overall gradient between generally hotter surface and cooler core of the body, but given the radiation source (Sun) has a specific direction with respect the body, this surface-core gradient must always be accompanied with the appropriate latitudinal and/or longitudinal temperature gradients.
... have
Rearranging the terms, we can express this condition by saying that the S- and P-wave velocities and for signal (seismic) propagation in the material are always much larger than the linear circumferential speed at equator . Indeed, and are of the order of kilometers per second, being of the order of meters per second.

## All Figures

 Figure 1: Griffith's criterion for the brittle fracture of a material with a tensile strength  MPa: and are the maximum and the minimum eigenvalues of stress tensor - . Note the positive values of  describe a compression, while the negative values describe a tension. The material (ordinary chondrite in this case) fails above, and to the left of the critical curve shown in the figure. Open with DEXTER In the text

 Figure 2:Thermal stress and temperature distribution in a meridional section of a centimeter-size CC meteoroid at 0.14 AU heliocentric distance. The Sun (indicated by the arrow) is located above the north pole; as a result, the time-dependent component of both fields vanishes. Isolines show given constant values of the ( left panel) and ( middle panel) principal components of the stress tensor, both in MPa, and the temperature in K ( right panel). Open with DEXTER In the text

 Figure 3:Thermal stress and temperature distribution in a meridional section of a centimeter-size CC meteoroid at 0.14 AU heliocentric distance. The Sun (indicated by the arrow) is located above the equator) and the rotation frequency is assumed to be 100 Hz. Isolines show given constant values of the (left panel) and (middle panel) principal components of the stress tensor, both in MPa, and the temperature in K (right panel). Open with DEXTER In the text

 Figure 4: Thermal stress and temperature distribution for a slowly rotating centimeter-size CC meteoroid at 0.14 AU heliocentric distance. The Sun (the direction of which is indicated by the arrow in the left panels) is located at the colatitude with respect the spin axis (oriented upwards) and the rotation frequency is assumed to be 0.1 Hz. Left panels show a meridional section of the body (containing the solar direction), right panels show an Aitoff projection of the surface. Isolines show given constant values of the (top panels) and (middle panels) principal components of the stress tensor, both in MPa, and the temperature in K (bottom panels). Open with DEXTER In the text

 Figure 5: The principal component - of the stress tensor at the center of a meteoroid as a function its size and the heliocentric distance. The isobars are labeled according to their values in MPa. The upper plots for the solar colatitude and bottom ones for . The left panels for the CC material parameters, the right panels for the OC material parameters. The rotation frequency f is inversely proportional to the size D as in Eq. (68). The shaded zone indicates where the critical line of the Griffith criterion has been overrun and the thermal fission occurs. Open with DEXTER In the text

 Figure 6: The principal component - of the stress tensor at the center of a meteoroid as a function its size and the solar colatitude . The isobars are labeled by values in MPa. The left panel for the CC material parameters, the right panel for the OC material parameters. The rotation frequency f is inversely proportional to the size D as in Eq. (68) and the heliocentric distance is 0.14 AU. The shaded zone indicates where the critical line of the Griffith criterion has been overrun and the thermal fission occurs. Open with DEXTER In the text

 Figure 7: The principal component of the stress tensor in the center of a meteoroid as a function its size and the rotation frequency f. The isobars are labeled by values in MPa. The left panel for the CC material parameters, the right panel for the OC material parameters. The heliocentric distance is 0.14 AU and the chosen value of the solar colatitude maximizes the thermal stresses (Fig. 6). The solid line shows the f(D) relation from Eq. (68). The dashed line indicates the situation in which the radius of the body is five times the penetration depth of the diurnal thermal wave (Eq. (69)). The shaded zone indicates where the critical line of the Griffith criterion has been overrun and the thermal fission occurs. Open with DEXTER In the text