Issue 
A&A
Volume 519, September 2010



Article Number  A75  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201014281  
Published online  15 September 2010 
Thermal stresses in small meteoroids^{}
D. Capek^{1}  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 metersize, 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 planetcrossing 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 planetcrossing region (e.g., meteoroids or small nearEarth 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 nonuniform temperature field in the solid body that, under some circumstances, may exceed the material strength. Either crack formation or a breakup 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 preexisting microcracks leading eventually to the breakup. 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 nonicy planetesimals residing on highly eccentric orbits can be destroyed by thermal stresses in the neighborhood of UX Orionistype young stars producing a heavyelementrich 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 builtup of a lowconductivity surface (regolithlike) 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 breakup 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 coldmaintained core and the solarheated 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 nonradial 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 quasistatic 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 diurnalcycle approximation from (i) has also a
consequence on the assumed maximum size of the meteoroid. This is
because the modeled insolation field has all meanmotioninduced
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)
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 < 14 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
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
where is the thermal emissivity, W m^{2} K^{4} is the StefanBoltzmann 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 fourthorder, nonlinear term in Eq. (3). The traditional wayaround 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 higherorder 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 righthand side of
Eq. (3) can be expressed using a coupled series of
spherical harmonics in parameters
and a Fourier
series in time t
where is the solar flux at the given heliocentric distance and is the colatitude of the solar direction in the bodyframe system introduced above^{}. We recall that the diurnal approximation adopted in this paper implies that is a constant, timeindependent parameter. The solar colatitudedependent functions b_{nk} are given by (e.g., Vokrouhlický 2006)
b_{00}  (5)  
b_{10}  (6)  
(7)  
b_{nk}  (8) 
where the last row applies to , and P_{n} and P_{n}^{k} are Legendre polynomials and the associated Legendre functions, respectively. We note that P_{n}(0)=0 for n odd such that only evendegree terms appear in the insolation (4) for . Thus, the dipole term with n=1 is the only imprint of the northsouth asymmetry.
Similarly to the insolation term, the temperature field can also
be expressed using a mixed Fourier and sphericalharmonics
expansion
We note we are mainly interested in the temperature variations because only they can result in thermal stresses. The amplitude functions t_{nk}(r,t) can again be given in Fourier series
where j_{n}(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 j_{n}functions for the timedependent terms. A simple analysis shows that the latter terms show a quasiexponential 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 C_{n0} and C_{nk} in
Eqs. (10) and (11) have to be determined from the
boundarycondition constraints
(3). One directly obtains
(see also Vokrouhlický 2006,1998,1999)
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 nonzonal 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
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 TT_{0} with respect some reference value T_{0}. Placing all the terms together, we have a generalized Hook's law in the form
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 T_{0}, 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 DuhamelNeumann relation (e.g., Turcotte & Schubert 2002; Kupradze et al. 1979; Boley & Wiener 1960)
which for vanishing volumic forces (such as the fields of solar tides or inertial forces in a rotating body) reduces to a simpler form
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)
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 DuhamelNeumann 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).
Because of the linearity of the DuhamelNeumann equation and the convenience with which we can develop the thermal gradients on its righthand 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 wellknown 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 righthand 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 lefthand side of the DuhamelNeumann 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 nonzonal modes represent timedependent (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)
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 stressfield 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.050.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 chondritelike (OC for short) material analogous to the group I meteoroids; and (ii) primitive carbonaceous chondritelike (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 roomtemperature measurements of ,
K, ,
,
and
of a H5class 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 chondritelike material). We ensured that
there is a good match between the roomtemperature 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 temperaturedependence measurements of the forsterite's specific heat capacity c (J K^{1}kg^{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^{1}kg^{1} for T=300 K, while Medvedev et al. (1985) obtained c=830 J K^{1}kg^{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:
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
temperaturedependence measurements, was approximated with the
temperature dependence for the shear strength
(Ohnaka 1995; Rocchi et al. 2004)
where Q is the activation energy, is the gas constant, and T_{1} is the reference temperature. We adopt K and T_{1}=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 Stype 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 observationallydetermined 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
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 nonaxisymmetric terms are timedependent and correspond to a frequency for the order k in the spherical harmonics functions.
3.1.1 Timeindependent 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
where x=r/R and
we suppressed the parametric dependence of the radial profile Gfunctions keeping note only of its rdependence. With these functions defined, we then have
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 x^{n2} decays for small r. We should also point out that because its amplitude G_{2} is proportional to , the quadrupole contribution vanishes near the node of the seconddegree 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 (1x^{2})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 Rdependence 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 timedependent 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 timedependent part of the field is proportional to the coefficients that vanish for and .
3.1.2 Timedependent part of the thermal stress tensor
We now consider the timedependent part of the deformation and the stress field, which is expressed by the nonzonal () 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
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 temperaturedependent 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 DuhamelNeumann equation, given in Sect. A.2. All of them share the sphericalharmonics development and Fourierdevelopment 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 Q^{1}_{nk} and Q^{2}_{nk} need to be determined to satisfy the surface boundary conditions (19) and Sect. A.3. We then spectrally decompose the stresstensor field
where
The explicitly temperaturedependent 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 nonzero 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 righthand 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 DuhamelNeumann equation.
For the first mode, we obtained
Figure 2: Thermal stress and temperature distribution in a meridional section of a centimetersize CC meteoroid at 0.14 AU heliocentric distance. The Sun (indicated by the arrow) is located above the north pole; as a result, the timedependent 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 
We next indicate the thermal stress components of the particular solution of the DuhamelNeumann equation: . We find it suitable to define^{}
(56) 
and the auxiliary radial functions
(57) 
With them, we have
Finally, we recall that the explicit temperature contribution in the generalized Hook's law is a pure volumic expansion
Equipped with these partial results, we may now compute the coefficients Q^{1}_{nk} and Q^{2}_{nk}. 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
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 (n1)=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 Q^{1}_{nk} and Q^{2}_{nk} constants. After solving the
corresponding set of algebraic equations, we obtain
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 thermalstress 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 centimetersize 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 timedependent 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 timedependent 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 solardirection 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 centimetersize 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 timedependent 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 stresstensor 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 Griffithcriterion line (Fig. 1).
Nevertheless, we are left with the fourdimensional 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
where D is in meters and f in Hertz. This relation of inverseproportionality matches very roughly suggested rotation periods of centimeter to decimetersize 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 OCtype 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 takeaway 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 nonzero 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 thermalstress 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 timedependent 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 timedependent part of the stress field is
confined to a tiny surface layer and the frequencyindependent
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 timedependent
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
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 timedependent 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 timevariability 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 YarkovskyO'KeefeRadzievskiiPaddack (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 selfregulating 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 preatmospheric 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.30.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 groupIIIb meteors of an estimated meterrange preatmospheric 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 preatmospheric 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 seasonalthermal 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 longperiod 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 largestobserved 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. Seminumerical 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.
AcknowledgementsWe 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
 Anderson, O. L., Isaak, D. L., & Oda, H. 1991, J. Geophys. Res., 96, 18037 [NASA ADS] [CrossRef] [Google Scholar]
 Beech, M., & Brown, P. 2000, Planet. Space Sci., 48, 925 [NASA ADS] [CrossRef] [Google Scholar]
 Beech, M., & Nikolova, S. 2001, Planet. Space Sci., 49, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Boley, B. A., & Wiener, J. H. 1960, Theory of thermal stresses (New York: Wiley) [Google Scholar]
 Borovicka, J. 2007, in IAU Symp., ed. A. Milani, G. B. Valsecchi, & D. Vokrouhlický, 236, 107 [Google Scholar]
 Britt, D. T. & Consolmagno, G. J. 2003, Meteoritics & Planetary Science, 38, 1161 [Google Scholar]
 Brown, P. G., Hildebrand, A. R., Zolensky, M. E., et al. 2000, Science, 290, 320 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bullen, K. E. 1975, The Earth's density (London: Chapman and Hall) [Google Scholar]
 Ceplecha, Z. 1996, A&A, 311, 329 [NASA ADS] [Google Scholar]
 Ceplecha, Z., Borovicka, J., Elford, W. G., et al. 1998, Space Sci. Rev., 84, 327 [NASA ADS] [CrossRef] [Google Scholar]
 Capek, D., & Vokrouhlický, D. 2004, Icarus, 172, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Capek, D., & Borovicka, J. 2009, Icarus, 202, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Flynn, G. J. 2004, Earth, Moon, and Planets, 95, 361 [Google Scholar]
 Gounelle, M., Spurný, P., & Bland, P. A. 2006, Meteoritics & Planetary Science, 41, 135 [Google Scholar]
 Grinin, V., Natta, A., & Tambovtseva, L. 1996, A&A, 313, 857 [NASA ADS] [Google Scholar]
 Hall, K. 1999, Geomorphology, 31, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Iesan, D., & Scalia, A. 1996, Thermoelastic Deformations (Dordrecht: Kluwer Academic Publishers) [Google Scholar]
 Jenniskens, P. 2006, Meteor Showers and their Parent Comets (Cambridge University Press, Cambridge) [Google Scholar]
 Kaula, W. M. 1968, An introduction to planetary physics  The terrestrial planets (New York: Wiley) [Google Scholar]
 Kuehrt, E. 1984, Icarus, 60, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Kupradze, V. D., Gegelia, T. G., Basheleishvili, M. O., & Burchuladze, T. V. 1979, Threedimensional problems of the mathematical theory of elasticity and thermoelasticity (New York: NorthHolland Pub. Co.) [Google Scholar]
 Marchi, S., Delbó, M., Morbidelli, A., Paolicchi, P., & Lazzarin, M. 2009, MNRAS, 400, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Medvedev, R. V., Gorbatsevich, F. I., & Zotkin, I. T. 1985, Meteoritika, 44, 105 [NASA ADS] [Google Scholar]
 Oberst, J., Molau, S., Heinlein, D., et al. 1998, Meteoritics & Planetary Science, 33, 49 [Google Scholar]
 Ohnaka, M. 1995, Geophys. Res. Lett., 22, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Opeil, C. P., Consolmagno, G. J., & Britt, D. T. 2010, Icarus, 208, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Parkus, H. 1976, Thermoelasticity (Wien  New York: SpringerVerlag) [Google Scholar]
 Paterson, M. S., & Wong, T. 2005, Experimental Rock Deformation  The Brittle Field (Berlin: Springer) [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes, 2nd edn (New York: Cambridge University Press) [Google Scholar]
 Rocchi, V., Sammonds, P. R., & Kilburn, C. R. J. 2004, Journal of Volcanology and Geothermal Research, 132, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Shestakova, L. I., & Tambovtseva, L. V. 1997, Earth, Moon, and Planets, 76, 19 [Google Scholar]
 Tambovtseva, L. V., & Shestakova, L. I. 1999, Planet. Space Sci., 47, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Tauber, F., & Kuhrt, E. 1987, Icarus, 69, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Thorne, K. S. 1980, Rev. Mod. Phys., 52, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Turcotte, D. L., & Schubert, G. 2002, Geodynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Urquhart, M. L., & Jakosky, B. M. 1997, J. Geophys. Res., 102, 10959 [NASA ADS] [CrossRef] [Google Scholar]
 Vokrouhlický, D. 1998, A&A, 335, 1093 [NASA ADS] [Google Scholar]
 Vokrouhlický, D. 1999, A&A, 344, 362 [NASA ADS] [Google Scholar]
 Vokrouhlický, D., & Bottke, Jr., W. F. 2001, A&A, 371, 350 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vokrouhlický, D. 2006, A&A, 459, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vokrouhlický, D., Nesvorný, D., Dones, L., & Bottke, W. F. 2007, A&A, 471, 717 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yomogida, K., & Matsui, T. 1983, J. Geophys. Res., 88, 9513 [Google Scholar]
Online Material
Appendix A: Solution of the DuhamelNeumann equation
We outline the main steps needed to solve the DuhamelNeumann 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 sphericalharmonicstype
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 purespin vector
harmonics, we have
with
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 radialprofile amplitudes U_{nk}(r) and V_{nk}(r). The second component in Eq. (A.1), , is the toroidal component.
The spheroidal character of our source (temperature) term the DuhamelNeumann equation implies two simplifications. First, the toroidal part of the displacement vector becomes negligible and we have W_{nk}=0. Second, we can restrict the summation over degrees n in Eq. (A.1) to the dipole and higherorder 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 nontrivial, purely radial profile^{} and the only viable free monopole term must have U_{00}=W_{00}=0to match the boundary conditions. Finally, we note that we also anticipated the Fourierdevelopment structure in Eq. (A.1) as it follows from the source ( development).
Substituting the spheroidalvector representation of into the DuhamelNeumann Eq. (18) we obtain the following system of equations for the radial profile of the amplitude functions U_{nk}(r) and V_{nk}(r) ():
and
Here we found it useful to separate the r and tdependences of the t_{nk}(r,t) amplitudes of the development in Eq. (9) and introduce pure radial parts T_{nk}(r) such that . These represent source terms in Eqs. (A.4) and (A.5).
While the solution of U_{nk}(r) and V_{nk}(r) is coupled by means of Eqs. (A.4) and (A.5), the fundamental implication of the DuhamelNeumann 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 U_{nk}(r) and V_{nk}(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
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 orderindexes n and k)
The partial derivatives of the spherical functions Y_{nk} 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 stresstensor field. Because of the explicit analytical solution for , this is achieved at no computational expense.
Because of the linearity of the DuhamelNeumann 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 DuhamelNeumann equation
Equations (A.4) and (A.5) with zero righthand sides represent the homogeneous DuhamelNeumann 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 DuhamelNeumann equation become degenerate in the k (order) index of the sphericalharmonics 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 powerlaw
solutions
U^{i}_{nk}=Q_{i} r^{mi} and
V^{i}_{nk} = r^{mi} with
,
realvalued exponents m_{i} and amplitudes Q_{i}.
After a straightforward algebra, we obtain
U^{1}_{nk}(r)  =  (A.15)  
V^{1}_{nk}(r)  =  r^{n+1},  (A.16) 
U^{2}_{nk}(r)  =  (A.17)  
U^{3}_{nk}(r)  =  (A.18)  
U^{4}_{nk}(r)  =  (A.19)  
V^{4}_{nk}(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 DuhamelNeumann equation
We next find a particular solution of the inhomogeneous DuhamelNeumann equation with the thermal source . We divide this task into a discussion of the stationary case (k=0) and timedependent case (). In both cases, we again use the approximation of neglecting the terms in Eqs. (A.4) and (A.5).
A.2.1 Timeindependent part
The stationary temperature field is given by
(Eqs. (9) and (10)) and thus
.
We again search the fundamental
system of solutions in a powerlaw form
and
with some realvalued 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 Timedependent part
The timedependent 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 DuhamelNeumann equation
Eq. (18), we obtain
where we have suitably assumed that the arbitrary constant on the righthand 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 sphericalharmonic and Fourier structure of the source term on the righthand side. Assuming thus a separation
we obtain
and
where is the elastic Pwave 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 amplitudefunctions of the spheroidalfield 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 DuhamelNeumann
equation is a linear combination of the freespheroidal modes
and
from Sect. A.1
and the particular mode
from Sect. A.2.
In the individual spherical harmonics modes, we have
,
where Q^{1}_{nk} and Q^{2}_{nk} 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
or again in the spherical harmonics modes
The truly active and independent conditions are and , from which the two constants Q^{1}_{nk} and Q^{2}_{nk} 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 timedependent part of the stress field, the algebra is more involved and we could not reach as simple and compact results as for the timeindependent part. We thus confine ourselves to provide formulae for the stresstensor components of the individual components and those for the integration constants Q^{1}_{nk} and Q^{2}_{nk}(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 rotationfrequencydependent 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, meanmotion 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 solarilluminated 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 surfacecore 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 Pwave 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 centimetersize CC meteoroid at 0.14 AU heliocentric distance. The Sun (indicated by the arrow) is located above the north pole; as a result, the timedependent 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 centimetersize 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 centimetersize 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 
Copyright ESO 2010
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.