Charles El Mir

K T Ramesh

Marco Delbo

The efficiency of thermal fatigue in regolith generation on small airless bodies

Regolith generation by thermal fatigue has been identified as a dominant mechanism for the breakdown of small (cm-sized) rocks on certain airless bodies. Simple numerical models for thermal fatigue seemed to indicate that this breakdown occurs faster in the larger decimeter-sized rocks, which is in contrast to the predictions of disruption models through successive micrometeorite impacts. The observation is justified by the existence of larger temperature gradient in bigger rocks, but it is not clear that this conclusion can be extrapolated or scaled to meter-sized boulders. Here we reveal a transition in the rock disaggregation rates by thermal fatigue when rock sizes rise above a critical length scale. A simple analytic model is formulated to predict the time to fracture of rocks on small airless bodies. We consider an uncoupled approach consisting of a one-dimensional thermal model, and a two-dimensional fracture model. The solution of the heat equation is used as input to the thermomechanical crack growth problem.

understanding could provide bounds on the survival rates of asteroidal rocks, and may help in coupling thermal fatigue with a mechanical disruption model to obtain a multi-mechanism view of regolith evolution in the solar system.

Introduction

Asteroids, comets, and small planets in the solar system are covered with a layer of loose and unconsolidated rocks called regolith. Surface imagery acquired by space probes has indicated that even the relatively small kilometersized asteroids carry a complex blanket of fine rocks smaller than a few centimeters. Images of asteroid ( 25143) Itokawa [START_REF] Murdoch | Asteroid Surface Geophysics[END_REF] and comet 67P/Churyumov-Gerasimenko [START_REF] El-Maarry | Fractures on 53 comet 67P/Churyumov-Gerasimenko observed by Rosetta/OSIRIS[END_REF] highlighted the existence of freshly exposed boulder surfaces, indicating that the evolution of these boulders is driven by an active surface rejuvenation process that could differ from those manifested on the lunar surface (Veverka et al., 2001b).

It is generally accepted that the lunar regolith is the byproduct of a long history of mechanical disruption wherein large boulders are eroded by a series of micro-and macro-meteoritic impacts [START_REF] Gault | Mixing of the lunar regolith[END_REF]. The abundance of impact-induced agglutinates points towards the dominant role that meteoritic impact had on the evolution of the lunar regolith. However, on small asteroids, the ejecta velocities from an impact event typically exceed the asteroid's escape velocity, and most of the resulting small debris would be lost into space as opposed to being reaccumulated on the surface [START_REF] Housen | Asteroidal regoliths[END_REF]. Early models of impact-induced regolith on airless bodies indicated that small asteroids (10 km or smaller) should only retain a negligible regolith layer [START_REF] Housen | Asteroidal regoliths[END_REF][START_REF] Housen | Regoliths on small bodies in the solar 55 system[END_REF][START_REF] Pettengill | Radar observations of asteroids[END_REF]. Consequently, asteroidal regolith generation remained, for the most part, a poorly understood mechanism [START_REF] Melosh | Planetary Surface Processes[END_REF].

Recently, thermal fatigue was suggested as a mechanism for in-place rock breakdown capable of generating fine-scaled regolith layers without any subsequent ejection. Observational evidence in favor of thermal fatigue in the solar system includes the work of [START_REF] Eppes | Cracks in Martian boulders exhibit preferred orientations that point to solarinduced thermal stress[END_REF] who collected orientation measurements from more than 1,800 cracks visible in nearly 1,500 rocks photographed by the Spirit rover during its journey on the Martian surface.

Their measurements showed that these cracks exhibit preferred orientations consistent with solar-induced thermal stresses. [START_REF] Dombard | Boulders and ponds on the Asteroid 433 Eros[END_REF] also showed that some boulders inside regolith "ponds" on 433 Eros have debris aprons, which were interpreted as the evidence that these rocks erode in place, likely due to thermal cracking. In another work, [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF] examined fragmentation induced by the diurnal temperature variations on asteroids. Their results indicated that thermal fatigue could play an important role in the generation of fine-grain regolith on small asteroids, and that the process is several orders of magnitudes faster than mechanical impact. [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF] also showed that thermal fatigue could break down larger rocks faster than smaller ones, owing to the greater temperature gradients in the bigger rocks-a trend that is in contrast to what is seen in the case of mechanical disruption by microimpacts. This latter result provides an explanation for the presence of fine regolith on small asteroids, but it is not clear how it would apply for larger boulders (the [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF] calculations were for surface rocks with diameters between 1 and 10 cm). The presence of meter-sized rocks on several asteroids, such as those pictured on Itokawa (Veverka et al., 2001b;[START_REF] Saito | Detailed Images of Asteroid 25143 Itokawa from Hayabusa[END_REF], indicates that thermal fatigue ought to behave differently at these larger sizes.

Here we demonstrate that thermal fatigue has an inherent lengthscale that characterizes the bounds on the rate of fragmentation for rocks in the range of 1-50 cm. We first develop a thermomechanical model that tracks the evolution of a crack due to the diurnal temperature variations in an asteroidal rock and efficiently bridges across the vastly different timescales of rotation and time to fracture. Using the model, we compute the trends in the timeto-fragmentation of sub-meter-sized rocks. We then develop a generalized analytical expression to predict the survival times of regolith particles on small airless bodies. Note that the regolith typically consists of rocks of a wide range of sizes. Here we focus on rocks that are larger than a mm. The rock breakdown study performed in this paper is limited to the aforementioned range of 1-50 cm, and it will be shown that thermal fatigue would be most efficient for rocks of such diameters.

Numerical Modeling of Thermal Fatigue

Thermal fatigue cracking of rocks results from the fluctuations in the temperature field that develop driving forces (which are not merely the stresses, but rather the stress intensity factors at the crack tips) on cracks within a rock. These cracks may then experience sub-critical fatigue growth until fragmentation. In this section, we present a numerical framework that calculates the crack tip driving forces due to thermal cycling and predicts crack growth over a large number of cycles.

Thermal Model

We consider an uncoupled solution that allows us to solve the heat equation (solar radiation and conduction) first, and use the resultant temperature field as an input to the thermomechanical deformation problem. Effectively, we are assuming that the stress state does not affect the thermal properties of the rocks and that the cracks have negligible effects on the temperature field. The motivation for this assumption is that heat is conducted parallel to the direction of the crack propagation. We recognize that the presence of cracks affects the thermal conductivity of a rock, essentially in the direction perpendicular to the crack plane.

The temperature profile in surface rocks is calculated using a diffusion model based on the asteroid thermal model developed by Spencer et al. [START_REF] Spencer | Systematic biases in radiometric diameter determinations[END_REF]. This model calculates the temperature T (t, z) as function of time and depth in a rock at the equator of an asteroid by solving the heat equation:

ρC p ∂T ∂t = κ ∂ 2 T ∂z 2 , ( 1
)
where κ is the thermal conductivity, C p the heat capacity, ρ the density and the coordinate z increases from zero at the surface downwards. Numeric values for the parameters used in this study are given in Table 1. It is documented that thermal conductivity varies with temperature, in particular in the case of meteorites [START_REF] Opeil | The thermal conductivity of meteorites: New measurements and analysis[END_REF][START_REF] Opeil | Stony meteorite thermal properties and their relationship with meteorite chemical and physical states[END_REF]. However, this effect is small for the diurnal temperature variations considered in this study. Recent astronomical observations showed this effect also in case of asteroids [START_REF] Rozitis | Observing the variation of asteroid thermal inertia with heliocentric distance[END_REF]; with some potential hints also observed on the asteroid 6 Hebe by [START_REF] Marsset | 3D shape of asteroid (6)˜Hebe from VLT/SPHERE imaging: Implications for the origin of ordinary H chondrites[END_REF]). This effect is probably due to the heat conduction in a regolith that is dominated by the radiative term, which is proportional to T 3 . For the case treated here concerning rocks where heat is essentially transported by the solid heat conduction, however, a temperature dependent thermal conductivity would only slightly enhance the maximum diurnal temperature variations by some degrees (as shown in [START_REF] Delbo | Asteroid Thermophysical Modeling[END_REF]). Therefore, in this formulation we assume that the thermal properties do not depend on depth and temperature. As a consequence, the thermal cracking effects hereafter estimated are somewhat conservative. Given the uncertainties in all other thermo-mechanical parameters, we prefer to use the simpler model of temperature independent conductivity, which gives results

that can be compared with other literature studies.

The heat diffusion equation is solved with the surface boundary condition of

(1 -A)S r 2 µ(t) = σ rad ir T 4 z=0 -κ ∂T ∂z z=0 (2)
where: σ rad is the Stefan-Boltzman constant; ir the infrared emissivity, which is assumed to be equal 0.9 as common in studies of asteroid surfaces (see [START_REF] Delbo | Asteroid Thermophysical Modeling[END_REF] and references therein); A is the bolometric Bond albedo, which is typically a few percent for the generally dark asteroid surfaces, r is the heliocentric distance of the asteroid in au; and S is the solar constant at 1 AU of 1370 W m -2 . The parameter µ is related to the direction cosine between the normal to the surface element and the direction to the sun θ, with µ(t) = cos [θ(t)] for -π/2 < θ < π/2 and µ(t) = 0 otherwise. The angle θ represents the local day time of the surface element, i.e., θ = 2π/P where P is the rotation period of the asteroid. This boundary condition implies that the rock is treated as a flat surface. It is known that surface roughness affects the observed temperatures of airless bodies [START_REF] Rozitis | Directional characteristics of thermal-infrared beaming from atmosphereless planetary surfaces -a new thermophysical model[END_REF][START_REF] El-Maarry | Fractures on 53 comet 67P/Churyumov-Gerasimenko observed by Rosetta/OSIRIS[END_REF]; however, this is a global effect on large surfaces. Locally, roughness can enhance the dayside temperature of rocks on the surface due to the mutual heating. This happens when a rock has a nonzero view factor towards warm surface elements (i.e.: the rock experiences local radiation due to the local topolgy), resulting in slightly larger temperature excursion between day and night. The flat surface thermal model we consider here is thus a conservative case.

At a depth of 5 meters, the additional boundary condition of ∂T ∂z z=5m = 0 implies that there is no heat flow towards (or from) deeper into the body.

These are similar to the conditions imposed by classical asteroid thermal modeling (see [START_REF] Delbo | Asteroid Thermophysical Modeling[END_REF] for a review).

The resulting time-dependent and depth-dependent temperature fields for rocks on the surface of a small carbonaceous chondritic near-Earth asteroid having a rotational period of 6 hours are shown in Fig. 1. This is used as the input to the computational model in Section 2.4.

In Fig. 1, the temperature as a function of depth during a 6 hour day/night cycle is plotted at 15 min increments (changing hues from red to blue).

The asteroid's surface experiences the largest temperature changes, which are quickly attenuated and disappear below a depth of 20 cm centimeters.

Notice that the temperature variations as a function of depth are highly time-dependent, with the maxima being experienced at different times for different depths. This highlights the need for an efficient thermomechanical model that can solve for the stress fields at small increments of times, while being capable of tracking the crack growth for many years.

Thermomechanical Crack Growth Model

In this section, we develop a thermomechanical model that predicts the fatigue crack growth in a rock as a result of the spatiotemporally varying temperature field calculated in the thermal model (Section 2.1). To do this, we must compute the driving force on a crack tip in the rock. This driving force is called the stress intensity factor or SIF. It is important to recognize that the stress field computed within an uncracked body does not identify whether or not a crack will grow [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF], and indeed, an accurate analysis of the stress field within an elastic rock containing a crack will always show an infinite stress at the crack tip. Understanding crack growth requires the use of the concepts of fracture mechanics. In general, cracks grow in one of three modes, each of which has an associated stress intensity factor: Mode I or tension, which results from tension perpendicular to the crack faces, with stress intensity factor K I ; Mode II or shear, which results from shear stresses acting parallel to the crack faces and causing crack sliding, with stress intensity factor K II ; and Mode III or antiplane shear, which results from shear displacements acting parallel to the crack front, with stress intensity factor K III . These stress intensity factors fully describe the stress, strain and displacement functions at the crack tip, and crack growth laws (crack direction, crack length and crack speed) are typically written in terms of these quantities (S. M. [START_REF] Beden | Review of Fatigue Crack Propagation Models for Metallic Components Review of Fatigue Crack Propagation Models for Metallic Components[END_REF][START_REF] Paris | Service load fatigue damage -a historical perspective[END_REF]. Having multiple modes active at one time is said to define a mixed-mode problem, and the thermal fatigue problem turns out to be mixed mode. So, to model the thermal fatigue fragmentation caused by the growth of a crack in a rock, we must first compute the stress intensity factors as the temperature field is varying.

Background

Stress intensity factors depend on both the loading configuration as well as the geometry and crack size [START_REF] Tada | The Stress Analysis of Cracks Handbook, Third Edition[END_REF]. Closed-form analytical solutions only exist for some simple idealized cases such as remote tension on a cracked plate. Some "universal" weight function approaches have also been developed [START_REF] Glinka | Universal features of weight functions for cracks in mode I[END_REF] but their application remains limited, especially with increasing complexity in loading profiles and crack configurations. Consequently, numerical approaches such as the Finite Element Method (or its extensions) are often needed. While traditional finite element method formulations are not capable of capturing the discontinuity in the displacement field that is caused by the presence of a crack, some measure of the singularity at the crack tip can still be achieved by explicitly meshing a sharp notch as part of the shape model. Post-processing techniques can then be used to estimate the stress intensity factor either by fitting the displacement or stress field solutions to the near-crack analytical solution, or by means of an energy integral approach [START_REF] Chan | On the finite element method in linear fracture mechanics[END_REF]. However, the post-processing approach requires an increasingly fine mesh near the notch to obtain satisfactory results, which increases the computational cost and data preparation effort for each simulation.

As a consequence, numerical approaches that extend the traditional finite element method have been developed for modeling cracked bodies. These are often achieved by formulating special types of elements that include a singu-larity in their shape functions [START_REF] Blackburn | Calculation of Stress Intensity Factors at Crack Tips Using Special Finite Elements[END_REF]Benzley, 1974;[START_REF] Nash Gifford | Stress intensity factors by enriched finite elements[END_REF] and using them near the notch to reduce the required amount of fine-meshing. Other methods include the quarter-point finite element [START_REF] Henshell | Crack tip finite elements are unnecessary[END_REF][START_REF] Barsoum | On the use of isoparametric finite elements in linear fracture mechanics[END_REF], the boundary collocation method [START_REF] Bowie | A modified mapping-collocation technique for accurate calculation of stress intensity factors[END_REF][START_REF] Newman | An Improved Method of Collocation for the Stress Analysis of Cracked Plates with Various Shaped Boundaries[END_REF], and the boundary elements method [START_REF] Alarcon | The boundary element method in elasticity[END_REF][START_REF] Cruse | Boundary Element Analysis in Computational Fracture Mechanics[END_REF]. These techniques can reasonably capture the singularity for a stationary crack tip. However, any crack propagation would then require the regeneration of a new mesh that conforms with the new crack geometry. This "remeshing" step is often the most computationally expensive step in an analysis, making these methods less suitable for crack evolution simulations.

One technique to avoid this remeshing step is the so-called eXtended Finite Element Method (XFEM) (Belytschko and Black, 1999;[START_REF] Karihaloo | Modelling of stationary and growing cracks in FE framework without remeshing: A state-of-the-art review[END_REF]. Belytschko and Black (1999) first introduced XFEM by discretizing a crack in the mesh through the enrichment (provision of additional degrees of freedom) of elements near the crack tip and crack faces.

The enrichments exploit the partition of unity property inherent in finite elements [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] to include the asymptotic displacement field resulting from the crack without modifying the existing mesh. [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] then introduced the generalized Heaviside step function to characterize the discontinuous field across the crack faces away from the crack tip, and [START_REF] Daux | Interpretation of thermal emission. I. The effect of roughness for spatially resolved atmosphereless bodies[END_REF] detailed a procedure for multiple branched cracks through the use of a junction function.

In XFEM, cracks can be incorporated into an existing finite element mesh by first identifying the elements that contain the crack, and labeling their respective nodes as face-enriched nodes related to elements that are split by a crack, or tip-enriched nodes related to elements containing a crack tip (see Fig. 2). Cracks can then move freely within the mesh without requiring a reconstruction of the mesh itself, but rather a simple identification of the elements to be "enriched" with additional degrees of freedom. The displacement field can then be computed using the standard finite element solvers [START_REF] Bathe | Representation of singularities with isoparametric finite elements[END_REF]. The stress intensity factor is then calculated in the post-processing stage using an interaction integral method [START_REF] Nash Gifford | Stress intensity factors by enriched finite elements[END_REF][START_REF] Wilson | The use of the J-integral in thermal stress crack problems[END_REF].

Model Details

Here, we develop a custom XFEM implementation (based on [START_REF] Bordas | An extended finite element library[END_REF]) that solves for the displacement field in cracked rocks subjected to cyclic thermal loading. Since the near-surface temperature variations change quickly with time (see Fig. 1), the model is optimized to efficiently bridge between the relatively small timesteps within a single period of rotation (using 15-minute increments) and the longer timescale of fragmentation (10 3 -10 6 years).

Consider a cracked rock that is thermally cycled. We idealize the rocks as plane strain cylindrical 2D bodies. First, a representative mesh of the geometry (circular in our case) is generated irrespective of the cracks. Then, cracks are overlaid on the mesh, and the elements that are split by a crack and those that contain crack tips are identified (Fig. 2). These elements are "enriched"

with the additional degrees of freedom, so that the total displacement field u at a point x is written in terms of the regular FEM displacement (without a crack) and an XFEM displacement (related to the crack enrichments):

Figure 2: Mesh of a circular rock containing a surface crack. The gray-shaded elements are directly affected by the crack and will include additional shape functions (degrees of freedom). The nodes identified with a red square belong to elements that are "split" by the crack and will contain step-function enriched degrees of freedom. The nodes labeled with a green circle belong to the element that contains the crack tip and will be enriched with asymptotic functions to capture the singularity at the crack tip. Each element contains a sub-scale distribution of heterogeneities (inset on the right), which are assumed to be circular and distributed in a regular array. These inclusions will contribute to the global stress tensor. The left inset presents the crack-tip coordinate system and the domains for the interaction contour integral used to calculate the stress intensity factor. On the mesh, these contours span 3-5 elements around the crack tip. Symbols are defined in the subsequent pages of the main text.

u(x) = i N i (x)û i u F EM + j N j (x)ψ(x)a j u XF EM (3)
where û is the nodal displacement from traditional finite element formulation,

N the finite element shape function, the a j 's are the additional degrees of freedom related to the enrichment ψ that is asymptotic for crack tip enrichment [START_REF] Williams | The Bending Stress Distribution at the Base of a Stationary Crack[END_REF] and a Heaviside step-function for crack face enrichment.

The temperature field induces thermal strains, which are captured through additional nodal forces given for a unit thickness as:

f th = A BC th dA = A BC (α∆T ) IdA (4)
with B being the displacement differentiation matrix, C the stiffness tensor, th the thermal part of the total strain, α the thermal expansion coefficients, and I the identity tensor.

Thereafter, the global stiffness matrix and the force vectors are constructed in the same manner as the traditional finite element method [START_REF] Reddy | An introduction to the finite element method[END_REF]. The nodal displacements are then obtained using an implicit solver, and the gradients of the displacement field gives the total strain tensor .

The Cauchy stress tensor (σ) is computed from the mechanical strain m , using the linear constitutive relation:

σ = C : m = C : -th (5)
Asteroidal rocks contain inclusions with varying thermomechanical properties. The mismatch in thermal expansion coefficients between inclusion and matrix generates additional internal stresses as the matrix constrains the grain's expansion. An analogous behavior is often observed in composite materials subjected to high temperature variations [START_REF] Evans | The role of inclusions in the fracture of ceramic materials[END_REF][START_REF] Biernacki | Experimental study of large scale model composite materials under thermal fatigue[END_REF][START_REF] Ajaja | Damage accumulation in a carbon fiber fabric reinforced cyanate ester composite subjected to mechanical loading and thermal cycling[END_REF]. Our model is also capable of treating inclusions as heterogeneities that are part of the material mesh, but such analysis will not be addressed in this paper. In this work, the inclusions are treated as sub-scale (Fig. 2) to our mesh resolution. The effective contribution of the thermal expansion mismatch stress caused by a volume fraction f of inclusions is accounted for using a representative volume element approach following a modified Eshelby formulation [START_REF] Hsueh | Residual thermal stresses in ceramic composites. Part I: With ellipsoidal inclusions[END_REF] as:

σ T M = 1 -2ν inc E inc + 1 1 -f 1 + ν m 2E m + f 1 -f 1 -2ν m E m -1 (α m -α inc )∆T I (6)
where E and ν are the Young's modulus and Poisson ratio, respectively. The subscript m refers to the matrix (or bulk) material properties, and inc refers to properties of the average inclusion (such as a chondrule).

Now that we have the displacements, strains, and stresses, we can compute the stress intensity factor at the crack tip. Typically, the energy release rate at the crack tip (which is related to the stress intensity factor) can be extracted by means of a J-integral approach [START_REF] Rice | A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks[END_REF]). The J-integral is a path-independent line integral that is equal to zero for a simply connected closed loop around a region containing no singularities. When the integral path contains a crack, the J-integral is equal to the energy release rate for the crack growth, which is related to the stress intensity factor at the crack tip.

Note that the J-integral alone does not provide enough information to extract the individual measures of K I and K II in a mixed-mode crack problem, since it depends on their combined effects:

J = (K 2 I + K 2 II ) E * (7)
where E * is equal to the Young's modulus at the crack tip E for plane stress, and E * = E/ (1 -ν 2 ) for plane strain, with ν the Poisson's ratio at the crack tip. Note that K III is not relevant for these 2D problems.

In this work, we use the interaction integral approach (Yau et al., 1980;Belytschko and Black, 1999), which decouples the individual contributions of K I and K II by superposing the solution of assumed "auxiliary" fields (note that the linearity of the elastic problem implies that the superposition of two fields at equilibrium should lead to another equilibrium state). Proper choices of these auxiliary states allow for decoupling the K I and K II contributions.

The J-integral for the superposed equilibrium state (superscript sup below)

would contain contributions from the auxiliary state (superscript aux below)

and the original state in XFEM (no superscript). That is:

J sup = 1 E * tip (K I + K aux I ) 2 + (K II + K aux II ) 2 = J + J aux + I sup (8)
The I sup term is an interaction term that contains the coupling between the XFEM solution and the auxiliary field solution. It can be shown that I sup is a line integral of the form:

I sup = lim Γ 1 →0 σ aux ij u j,1 + σ ij u aux j,1 -σ aux jk jk δ 1i m i qdΓ 1 (9)
where m i are the components of the unit outward normal vector to the integrating contour Γ 1 (Fig. 2).

Recall that the fundamental property of the J-integral is that it equates to zero over a closed path away from the crack. However, [START_REF] Wilson | The use of the J-integral in thermal stress crack problems[END_REF] showed that, in general, the integral is non-zero when thermal stresses are present, meaning that for thermal stress crack problems the crack tip stress intensity factors cannot be determined directly from such a line integral calculation. The alternate formulation of [START_REF] Wilson | The use of the J-integral in thermal stress crack problems[END_REF] is used in our model after converting the line-integral into an equivalent domain integral (EDI) written in crack-tip coordinates [START_REF] Merle | Solving thermal and phase change problems with the eXtended finite element method[END_REF][START_REF] Kc | Interaction integrals for thermal fracture of functionally graded materials[END_REF], as this is better suited for implementation in finite element codes since it is written in terms of quantities readily available from the analysis step:

I sup = A σ aux ij u j,1 + σ ij u aux j,1 -σ jk aux jk δ 1i q ,i dA + A σ aux ij α(∆T ) ,1 δ ij qdA (10)
Here, q is a weight function with values varying smoothly from 1 on Γ 1 to 0 on Γ 0 (Fig. 2). The second integral comes from the thermal contribution to the J-integral. The different stress intensity modes can now be individually extracted through the proper choice of auxiliary fields (Yau et al., 1980). It follows from Eq. ( 8) that the mixed-mode stress intensity factors are related to the interaction integral as:

I sup = 2 E * (K I K aux I + K II K aux II ) (11)
Then, by equating Eqs. ( 11) and ( 10) and using the fundamental solution [START_REF] Sun | The Elastic Stress Field around a Crack Tip[END_REF] for a crack in pure mode-I (K aux I = 1, K aux II = 0), the K II terms are dropped and we are left with K I as the only unknown. Similarly, setting the auxiliary field to the pure mode-II solution (K aux

I = 0, K aux II = 1) gives K II .
Now that the mixed-mode stress intensity factors are known, we use an appropriate fatigue crack growth law to describe the rate of crack propagation after a thermal excursion cycle. We here use the so-called Paris' Law [START_REF] Paris | Service load fatigue damage -a historical perspective[END_REF], which relates the crack length, a, in a given number of cycles, N , to the excursion in the effective stress intensity factor, ∆K ef f :

da dN = C (∆K ef f ) n (12)
where C and n are material parameters fit to fatigue experiments. Note that these parameters generally depend on the material, environment, frequency, temperature, and loading. We use the energy release rate model [START_REF] Rhee | Mixed-mode stress intensity factor solutions of a warped surface flaw by three-dimensional finite element analysis[END_REF] to define the effective stress intensity factor, combining the K I and K II contributions:

K ef f = K 2 I + K 2 II ( 13
)
The crack is allowed to grow arbitrarily in the mesh, with an orientation determined by the maximum tangential stress criterion (Erdogan and Sih, With these equations implemented into our model, we are able to perform high-resolution simulations of thermal fatigue crack growth over multiple cycles (see Fig. 3). We further optimize our algorithm to make use of multiple processors, making it suitable to be run on high-performance computational clusters. The computational "bottlenecks" of the XFEM algorithm in this case arise from: a) assembling the global stiffness matrix and force vector, b) solving the set of coupled linear equations, and c) detecting the nodes for enrichment after cracks propagate. Therefore, we first parallelize the stiffness matrix and force vector assembly process through an algorithm that maps the element-wise contributions to their position in the global matrices. This allows for resolving elements from different parts of the mesh onto separate computational cores and populating their contribution to the global stiffness matrix and force vector using the map. In addition, the temperature variation in our thermal fatigue problem has no effect on the global stiffness matrix during a single cycle (when the current crack has not yet propagated), and only affects the global force vector (following Eq. ( 4)). So during a single step, we only build the global stiffness matrix once and send the force vector for each time increment onto a separate node for parallel computation (that is, all time increments can be solved simultaneously on different computational nodes; Fig. 3) using an efficient parallel implicit matrix solver [START_REF] Witkowski | Software concepts and numerical algorithms for a scalable adaptive parallel finite element method[END_REF]. When the crack propagates at the end of a thermal cycle, we detect affected elements using a sub-grid search algorithm that starts from elements surrounding the crack tip (instead of looping over all elements).

θ = 2 arctan   1 4 K I K II ± 1 4 K I K II 2 + 8   (14) 1
This typically reduces the search algorithm to a much smaller loop since the fatigue crack growth distance in a single cycle rarely exceeds a few elements from the original position. These features allow us to efficiently bridge across a large number of cycles and makes our thermomechanical model effective for modeling the thermal fatigue breakdown of rocks on airless bodies in the solar system.

Stresses in cracked rocks under thermal cycling

We first discuss the stress fields that arise in a cracked rock as a result of thermal cycling. As a typical example, Fig. 4 The stress field perpendicular to the crack faces (stresses in the crack opening direction) is shown on the second and third rows. The second row is for a rock without any crack, while the third row shows the stress field in a rock of the same size that has a 15 cm surface crack. Notice that the stress at the crack tip is singular and that the crack relaxes the stress field at its faces.

Rate of rock breakdown by thermal fatigue

We now use the thermal and thermomechanical models presented so far to determine the time required to fracture rocks of varying diameters through the growth of a crack by thermal fatigue. We consider the case of a near-Earth asteroid with a 6-hour period of rotation, which is typical of many small asteroids [START_REF] Gregersen | The Outer Solar System: Jupiter, Saturn, Uranus, Neptune, and the Dwarf Planets[END_REF]. We model circular rocks of different diameters, all containing an initial surface crack (see the mesh in Fig. 2 and Fig. 8 in [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF]). We use the thermomechanical properties for carbonaceous chondrite (Table 1) and fatigue parameters of Carrara marble [START_REF] Migliazza | Experimental investigation on crack propagation in Carrara marble subjected to cyclic loads[END_REF]. These are the same material properties and initial conditions as in [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF]. To the best of our knowledge, the only performed also an uncertainty analysis and showed that values different than more than 30% from those of Carrara marble would not explain the laboratory experiment data on meteorites. The stress intensity factor is recorded at time increments of 15 minutes. At the end of a complete cycle, the excursion between the maximum and the minimum stress intensity factor is determined and the equivalent stress intensity factor excursion is calculated using Eq. ( 13). Then, Eqs. ( 12) and ( 14) are used to calculate the crack length increment and direction.

W m -1 K -1 M L T -3 θ -1 0.5 1.88 Opeil et al. (2010) Specific Heat Capacity Cp J kg -1 K -1 L 2 T -2 θ -1 500 550 Opeil et al. (2010) Bulk density ρ kg m -3 M L -3 1662 3150 Opeil et al. (2010) Thermal inertia Γ J m -2 s -0.5 K -1 M T -5/2 θ -1 640 1800 Opeil et al. (2010) Rotation period τ h T 6 6 -
Although our modeling approach is able to solve for arbitrarily growing cracks, recall that the temperature field computed from the thermal model is one-dimensional (that is, the temperature at a given time is a function of only one direction) and that the geometry and crack configuration are symmetric. This symmetry in geometry and loading implies that no mode-II fragmentation will occur (∆K II = 0 in every cycle), and consequently the crack will propagate over a straight path without any kinking (θ = 0). Given these idealizations, we consider the rock to be broken when the initial surface crack length becomes nearly equal to the rock diameter; that is, when the crack extends to the bottom-most element in the mesh. Note that a break in symmetry could occur when using a 2-D temperature field, or including multiple initial cracks, or an initially inclined surface crack. In these cases, fragments could also be produced by the coalescence of multiple cracks or by flaking surface material. While such scenarios can be captured reasonably through our numerical model, they are not considered in this work in the interest of first extracting some key physical ideas that can be simply applied for different asteroids in the solar system.

We have performed thermal fatigue simulations for such carbonaceous chondrites of diameters between 1 cm and 50 cm, and our results on the thermal fragmentation of these rocks are presented in 2014), but our more accurate results are still within their reported error margins. Our simulations also show that small rocks will take longer to fracture through thermal fatigue than larger rocks (up to approximately 7 cm). The trend captured here and in [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF] has also been stipulated by [START_REF] Molaro | Thermally induced stresses in boulders on airless body surfaces, and implications for rock breakdown[END_REF] for lunar rocks (albeit from purely stress-based arguments). However, Fig. 5 shows that beyond approximately 7 cm, this trend is reversed, in that larger rocks require increasingly longer times to fragment. This is primarily because although larger thermal gradients may cause larger stresses, the largest temperature fluctuations are limited to the near-surface (Figs. 1 and4), and temperature changes become negligible far from the skin depth. When cracks reach a size larger than the skin depth, the driving force at the crack tip gradually decreases, and the rate of growth of the crack is slowed down.

For comparison, Fig. 5 also shows the survival times of the rocks as a consequence of the micrometeorite impact mechanism, as calculated using the results from [START_REF] Hörz | Catastrophic rupture of lunar rocks: A Monte Carlo simulation[END_REF], with the modifications applicable to asteroid surfaces made in [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF]. In [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF], the mechanical disruption estimates were calculated for the 99% probability of survival.

However, [START_REF] Hörz | Catastrophic rupture of lunar rocks: A Monte Carlo simulation[END_REF] notes that the maximum survival times could be up to 25% shorter than this upper bound, depending on the exact particle abrasion rates and the continuously decreasing effective cross-section caused by the abrasion. To provide a more conservative comparison between thermal fatigue and mechanical disruption, we consider here the 50% probability of survival of the rocks by mechanical impacts (dashed line in Fig. 5). Note that observational evidence suggests that the rock breakdown rate on the Moon is at least a factor of 5 higher than the Horz model predicts [START_REF] Basilevsky | Survival times of meter-sized boulders on the surface of the moon[END_REF]. The model may thus underestimate breakdown rates due to micrometeorite impacts. Comparison of the thermal fatigue and impact breakdown mechanisms thus shows that the impact process may become at least as efficient a breakdown process as thermal fatigue for rocks as large as 50 cm, and is likely to be more efficient than thermal fatigue for larger rocks. Note again that this computation is for C-type asteroids at 1 AU and with a diurnal cycle of 6 hours.

The driving force for such sub-critical fatigue crack growth is the excursion in stress intensity factor ∆K, as shown in the Paris law (Eq. ( 12)). For However, as rock diameters become larger than the skin depth, a reversal in the trend is observed. That is, the crack growth speed of thermal fatigue driven surface crack is greatly reduced in these larger rocks.

all crack sizes shown here, the maximum stress intensity factor is always smaller than the fracture toughness K IC , and the crack grows only in fatigue and no dynamic fracture occurs.

In order to obtain a better understanding of the non-monotonic effect of thermal fatigue with respect to rock size, Fig. 6 shows the computed excursion in stress intensity factor (normalized by the material's fracture toughness)

as a function of crack length (normalized by the rock diameter) for rocks of a range of sizes. Note that when the normalized crack length reaches 1, the rock is completely broken. Looking first at the rocks that are less than 10 cm in diameter, we see that there is an increasing trend in ∆K with increasing crack length. This means that for these rocks, the driving force for crack growth increases as the crack grows (so that for the same thermal cycles, the rock is increasingly likely to break as the initial crack increases in length as a result of prior thermal cycles). Thus the crack growth accelerates as the crack grows. However, this is only true for rocks that are smaller than the diurnal skin depth. For the larger rocks, with diameters > 10 cm in this case, we see from Fig. 6 that the crack tip driving force ∆K decreases as the crack grows, leading to a decrease in the crack growth rate, and thus these larger rocks take an increasingly long time to break. This is true even though these larger rocks do sustain large thermal gradients. The decrease in the driving force is greater for larger rocks in this size domain, and thus larger rocks take a longer time to fragment through thermal fatigue. Essentially our detailed simulations show that for rock sizes significantly larger than the diurnal skin depth, thermal fatigue can efficiently grow cracks up to a certain size (roughly 60% of the diameter), but after that crack growth is greatly slowed down.

Note that this observation is largely independent of the chosen crack growth law. The crack tip driving force that we compute here in Fig. 6 depends on the current crack length rather than the crack growth history, and so this change in crack tip driving force will be observed regardless of the specific growth law. Here we assume a simple Paris-type relationship for fatigue growth, which may not properly describe the crack tip speed in thermal fatigue. Further experimental work is required to determine the applicability and the parameters of the Paris' law for thermal fatigue crack growth in different meteoritic materials. For a different growth law, the computed crack tip driving force will essentially be the same, but the computed rate of growth of the crack would be different (i.e., how long it takes the crack to cross the rock will change, but the trends with respect to the size effects will be the same).

The results of Fig. 6 indicate that there is a characteristic lengthscale that separates two distinct regimes in thermal fatigue: one that is dominated by larger temperature fluctuations happening near the crack tip, and another where the temperature changes become far-field relative to the crack tip location. This characteristic lengthscale is, to first order, determined by the diurnal skin depth (Fig. 1). In addition, the two trends described here and shown in Fig. 5 both appear to follow a power law-type dependence. That is, they may be fit to two distinct slopes on a log-log curve (in the same manner as the mechanical impact scaling). In the next section, we will use the concepts of this lengthscale and size-dependence of the excursion in stress intensity factor to derive a simplified relationship for the survival rate of rocks as a function of rock diameter.

Scaling Analysis

Our fracture mechanics model computes the stress intensity factor excursion during each asteroid rotation until the fragmentation of the rock. The rock's time to fracture is determined essentially by the number of thermal cycles that occur before the crack length becomes equal to the rock diameter. However, the amount of crack growth in a given cycle is not constant, so that the crack could (for example) spend a very long time (a large number of cycles) growing very slowly, and then accelerate rapidly and cover the remaining length in a very few cycles. Thus, in order to obtain a reasonable analytical estimate of the time it takes to crack a given rock, we primarily need to know the conditions under which the crack spends the largest fraction of its time, when it is growing very slowly. Now, the fatigue crack propagation law adopted in this work is the classical Paris law, relating the excursion in stress intensity factor at the crack tip to an incremental displacement of the crack front. In principle, it is therefore sufficient to integrate the Paris law to obtain an estimate for the predicted time to fracture. As Fig. 6 shows, the excursion in the stress intensity factor is itself a fairly nonlinear function of crack length, especially in the largest rocks. However, in terms of estimating the time to fracture, all we need is to know how ∆K depends on crack length during the times when the crack is growing very slowly, since this dominates the total time to fracture.

Our simulations show that for "small rocks," i.e., rocks smaller than the skin depth, ∆K increases approximately linearly with crack length for at least 80% of the total time before fracture (this occurs until the crack length reaches about 70% of the rock's diameter). The excursion in stress intensity for rocks smaller than the skin depth with normalized crack lengths a/d between 0.1 and 0.7 can be fit to a linear function ∆K(a) ≈ Āa + B with a Pearson correlation coefficient (Lee [START_REF] Rodgers | Thirteen Ways to Look at the Correlation Coefficient[END_REF]) r = 0.988. This linear approximation greatly simplifies the integration of the Paris law, Eq. ( 12) . In Section 4, we will discuss some of the consequences of this approximation. Now, using the assumption that ∆K(a) ≈ Āa + B, with Ā and B the linearization parameters, we can integrate the Paris fatigue crack growth law (Eq. ( 12)) from the initial crack length a 0 at time t = 0 to a final crack length a f = d at the approximated failure time t f . Now, recognizing that dN = 1/τ dt where τ is the period of rotation, we get:

t f τ ∼ 1 C Ā(1 -n) • Ād + B 1-n (15)
Note that here we only seek to understand the power-law dependence of t f on d. We do not calculate the exact value of the integral since the linearization parameters Ā and B are different for different rock diameters.

After some analytical simplifications (dropping the B terms and setting and setting t = t f at a = d) to Eq. ( 15), we can express the time to fracture as a function of rock diameter as:

d ∼ t f τ 1 1-n Ān C (n -1) 1 1-n (16a) t f τ ∼ (d) 1-n 1 Ān C(n -1) (16b)
Eq. (16a) and Eq. (16b) give us relatively simple approximate relationships between the time to fracture and the size of the rock (for rocks smaller than the diurnal skin depth), and shows that this relationship is dominated by the power 1 1-n in the former case and 1 -n in the latter case. Notice that t f /τ is dimensionless, and A n C has base units of L 1-n , making the right-hand side of Eq. (16b) also dimensionless.

For the case of rock diameters larger than the diurnal skin depth, the crack grows quickly to a normalized length of a/d = 0.2 in less than 10% of the total time to fracture. The bulk of the time to fracture occurs as the crack is growing progressively slower until a/d ∼ 0.95. For this domain, a linear fit does not capture adequately the shape of the excursion in stress intensity factors, and so the integration cannot be simplified in the same manner as for the smaller rocks. Instead, we directly fit to the final fragmentation times in Fig. 5 and obtain that d ∝ t 1 n-1 . We will show that these two slopes provide a good description of the thermal fatigue lifetime of rocks on asteroids with different composition, heliocentric distances, and rotation periods.

In order to understand the contributions of the thermophysical and mechanical variables to the final fragmentation time, we make use of the Buckingham π theorem [START_REF] Buckingham | On physically similar systems; Illustrations of the use of dimensional equations[END_REF]. The different π-groups are identified by first constructing the dimensional matrix [START_REF] Brand | The Pi theorem of dimensional analysis[END_REF] of all variables relevant to the problem (Table 1). Mathematically, the π-groups represent the null-space of the dimensional matrix that are obtained from the matrix's kernel vector. In this analysis, we have the following π-groups:

π 1 = d 2 • ρ • C p κ • τ (17) π 2 = d 4 • ρ • α κ • τ 3 (18) π 3 = κ • τ 3 • ∆T d 4 • ρ (19) π 4 = τ 2 • E d 2 • ρ (20) π 5 = t f τ (21) (22)
One must recognize, of course, that this set is not unique as any combination of the dimensionless groups also yields a dimensionless group. However, this particular set has some useful attributes. π 1 is effectively a measure of the diameter normalized by the skin depth, δ. Eq. ( 16b) can be re-expressed in dimensionless form after normalizing the diameters by the skin depth as:

d δ = Λ t f τ m ;          m = 1 1-n , d δ ≤ 1 m = 1 n-1 , d δ > 1 ( 23
)
where Λ is a scaling parameter related to the number of cycles needed to fragment a rock with diameter equal to the skin depth. In this form, the fragmentation time is expressed as a simple a power-law type relation with the rock diameter and skin depth. The advantage of this scaling relation is that the survival time of rocks by thermal fatigue can be estimated by performing a single XFEM simulation for a rock of size equal to the skin depth to determine Λ, and extrapolating to rocks with different sizes using Eq. ( 23).

This scaling analysis thus provides us with a simple way to understand the effects of thermal fatigue in terms of the time to fracture and rock size.

Using this simple scaling, Fig. 7 shows a normalized representation of Eq. ( 23) and the survival rates of rocks on small asteroids with the same rotational period but with different thermomechanical properties (both ordinary and carbonaceous chondrites) and heliocentric distances (both near-Earth asteroids and main-belt asteroids). In addition, we plot a normalized version of the data presented in [START_REF] Ravaji | The effect of rotation period on thermal stress weathering[END_REF]. Ravaji et al. calculated the predicted thermal fatigue lifetime of a 10 cm diameter ordinary chondrite-like rock on a near-Earth asteroid, as a function of the asteroid's rotation period using the same model of [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF]. They indicate that extrapolation of rock lifetimes for different rotational periods is not linear, but did not indicate how the scaling should be performed. In Fig. 7, we show that by scaling the rock diameter by the skin depth, which itself is a function of rotational period, the survival times of [START_REF] Ravaji | The effect of rotation period on thermal stress weathering[END_REF] collapse well onto the slopes identified in this work. Despite using a different thermomechanical model, these estimates fit well onto our nondimensional plot, especially for those rocks that are smaller than ∼ 4 skin depths. Note that beyond around 5 skin depths, the efficiency of thermal fatigue is greatly reduced, and the scaling arguments overestimate the breakdown rates by an order of magnitude.

What remains is to obtain an understanding of the meaning of the scaling parameter Λ. We do this by using a parametric analysis with our simulations to extract the dependence of Λ on C, E, α, and ∆T . This allows us to rewrite Eq. ( 23) as :

d δ ∼ C 1 1-n E -1.25 ρ -0.5 C -0.5 p (α∆T ) -1.5 t τ m (24)
This is a powerful relationship, because we can now estimate the survival rates of rocks for the same material and similar rotational periods with different positions in the solar system (that is, with different ∆T ).

As an example, by scaling the results from Figs. 5 and 7 using the scaling factors in Eq. ( 24), we present in Fig. 8 the predicted time to break down (defined by the colored contours) a 10 cm rock on a C-type asteroid as a function of the period of rotation and the heliocentric distance. The dark blue contours are of the order of 10,000 years, and the yellow corresponds to a million years (as shown in the key on the right of the figure). Note that the end of 2018. We hope that our results will therefore help in interpreting the observations from sample characterization for both missions.

We also note the existence of an intermediate range in rotational periods (4 -15 h) where our model predicts that thermal fatigue would be most efficient. Very fast rotators have small skin depths, and larger rocks would require increasingly larger times to fragment. Slow rotators have larger skin depths and rocks much smaller than that skin depth would also take increasingly longer times to fragment (see discussion for Fig. 5), as well as requiring a longer time to complete a single period. 24). Some C-type asteroids are highlighted on the plot using their semi-major axis and rotation period as reference. These were not modeled directly, and discrepancies from Table 1 in their particular surface and material composition would lead to different breakdown rates than those predicted in this figure. These results indicate the existence of an "intermediate" range in rotational periods where thermal fatigue would be more efficient. Very fast rotators have small skin depths, and larger rocks would require increasingly larger times to fragment. Slow rotators have larger skin depths and take longer time to complete a single period, and rocks much smaller than that skin depth would also take increasingly longer times to fragment (see discussion for Fig. 5). Consequently, the fragmentation times and the ideal "intermediate range" in rotation period depend on the particular rock diameter and how it compares with the skin depth.

Discussion

For cm-sized rocks, we predict that thermal fatigue can be several orders of magnitude faster than micrometeorite impacts in fragmenting such small rocks on small airless bodies with a ∼ 6 h period of rotation, in-line with the results of [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF]. In [START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF], thermal fatigue appeared to occur faster in larger rocks, owing to the larger temperature gradients they experience. Their results for ≤ 10 cm rocks indicated a speedup of thermal fatigue as rock sizes increase, implying that much larger rocks should be relatively scarce on the surface of asteroids. Similarly, Molaro et al.

(2017) reported an increase in peak stresses for large boulders (≥ 1 m) on the lunar surface as opposed to ≤ 30 cm rocks. They concluded, using purely stress-based arguments, that this could suggest a lack of very large boulders (diameters of several skin depths) on airless bodies. However, we have shown here that thermal fatigue is characterized by a signature lengthscale related to the thermal skin depth of the asteroid. While peak stresses due to thermal gradients can increase in larger rocks, they are mostly limited to the nearsurface and can be efficiently relaxed with the introduction of a crack. Peak stresses derived from purely elastic solutions are therefore poor indicators of the fatigue behavior. As cracks grow to lengths much larger than the skin depth, the driving force is greatly reduced, leading to a slow-down in thermal fatigue crack advance. Hence, our results imply the existence of a size domain over which thermal fatigue is mostly active.

These newfound size domain for thermal fatigue could help in determining the abundance and distribution of cm-sized rocks, which will be the main excavation targets in asteroid mining and sampling missions. For instance, it has been shown that asteroid (433) Eros has a complex regolith (Veverka et al., 2001b,a) whose size-frequency distributions [START_REF] Thomas | Shoemaker crater as the source of most ejecta blocks on the asteroid 433 Eros[END_REF] suggested that the mechanisms for depleting 20 cm diameter rocks seemed to be very different than those for 100 m blocks. Their observation is in-line with the lengthscales determined in this study, which may imply that thermal fatigue could be a contributing mechanism for the disaggregation of these relatively small surface rocks. In addition, the recently launched OSIRIS-REx mission that is targeting the asteroid (101955) Bennu [START_REF] Lauretta | OSIRIS-REx: Sample Return from Asteroid (101955) Bennu[END_REF] will attempt to sample and return a minimum of 60 g from the asteroid's bulk regolith. Bennu is a carbonaceous asteroid with a 4.3 h rotation period and an estimated thermal skin depth of ∼ 5 cm. In this work, we predict that for such a case, most cm-sized rocks should have experienced some degree of thermal fatigue, and we expect that the sampled rocks returned by OSIRIS-REx should show evidence of an active thermal fatigue mechanism.

During the preliminary survey of Bennu, in December 2018, OSIRIS-REx returned images covering 80% of the asteroid surface with a spatial resolution of 33 cm/px, visible, near and thermal infrared spectra [START_REF] Lauretta | The unexpected surface of asteroid (101955) bennu[END_REF]. Images show a geomorpholically diverse surface with craters (likely due to impacts) and covered with boulders of sizes ranging from 58 m of length down to the spatial resolution limit [START_REF] Della Giustina | The regolith of (101955) bennu from osiris-rex imaging and thermal analysis[END_REF][START_REF] Walsh | Craters, boulders and regolith of (101955) bennu indicative of an old and dynamic surface[END_REF]. Some of the boulders present evidence of fractures [START_REF] Walsh | Craters, boulders and regolith of (101955) bennu indicative of an old and dynamic surface[END_REF], while other boulder arrangement suggest that these are breaking in place [START_REF] Lauretta | The unexpected surface of asteroid (101955) bennu[END_REF]. Thermal fatigue is a potential explanation for the observed features. Another explanation is thermal dehyadration or desiccation of Bennu's boulders, or a combination of the two. OSIRIS-REx near-infrared spectra clearly indicate the presence of a 20-25% deep absorption band at 2.72µm, that is typical of hydrated meteorites of type CI or CM [START_REF] Hamilton | Evidence for widespread hydrated minerals on asteroid (101955) bennu[END_REF]. As such, it is possible that dehydration of this material could be responsible for the presence of fractures on boulders. The cumulative size distribution of Bennu's boulders larger than 8 m (which is the completeness limit) is well represented by a power law with index 2.9 ± 0.3. On the other hand, asteroid (25143) Itokawa has a value of 3.5 ± 0.1 for the exponent of the power law that best fits its cumulative size distribution of boulders ≥ 10 m [START_REF] Mazrouei | Block distributions on itokawa[END_REF]. The area number density of large boulders (≥ 20 m) on Bennu and Itokawa is comparable [START_REF] Della Giustina | The regolith of (101955) bennu from osiris-rex imaging and thermal analysis[END_REF] and references therein). This implies that more surface of Itokawa is covered in small particles than is the case for Bennu. In addition, OSIRIS-REx imaging data from Bennu's preliminary survey rule out a surface covered with centimeter or decimeter sized particles, but are consistent with the possible presence of much finer particles [START_REF] Della Giustina | The regolith of (101955) bennu from osiris-rex imaging and thermal analysis[END_REF]. In summary, OSIRIS-REx preliminary survey results supports fracturing process at the surface of Bennu and the scarcity of centimeter or decimeter sized particles, consistent with theoretical prediction from this work. At the beginning of March 2019, OSIRIS-REx will commence its detailed survey of Bennu, during which a global coverage of the surface at the spatial resultion of few centimeters will be obtained. The detailed survey data will enable a better understanding of the smaller rocks on Bennu to be obtained.

Furthermore, the effects of thermal fatigue could possibly be identifiable in the thermal inertia measurements. It has been shown that the thermal inertia of surfaces covered with small particles (sand-like regolith) is smaller than those covered by rocky fragments (of ∼ 10 cm diameter) [START_REF] Golombek | Size-frequency distributions of rocks on Mars and Earth analog sites: Implications for future landed missions[END_REF]. In that sense, a small fast-rotating asteroid with a large surface temperature variation would generate an increased number of small fragments, which could lead to a decrease in the measured thermal inertia.

Therefore, in the cases where thermal fatigue is driving the regolith evolution, thermal inertia (and skin depth) could be direct indicators of the evolution history of the rock sizes that blanket an asteroid's surface. Further research is needed to explore the possible relationship, if any, between thermal inertia and thermal fatigue.

For sufficiently large rocks, our model suggests that thermal fatigue would be more efficient in flaking off material from the surfaces as opposed to a complete fracture, which results in an apron of small (cm-sized) rocks surrounding larger boulders. This behavior has been observed on (433) Eros [START_REF] Robinson | The nature of ponded deposits on Eros[END_REF], where the aprons provided evidence that boulders are eroding in place, giving a probable explanation for the nearly 300 pond-like features on Eros [START_REF] Dombard | Boulders and ponds on the Asteroid 433 Eros[END_REF]. Such cases may be investigated through our model by solving for the 2D temperature field and tracking the growth of multiple crack networks in the thermomechanical XFEM model. The results derived here rely on a 1D thermal solution of an idealized spherical asteroid.

The temperature profile does not treat surface roughness, shadowing effects, or changes in thermal conduction introduced by cracks or rock boundaries.

A more sophisticated thermal model may lead to a different temperature profile than that in Fig. 1, which would give a different crack tip driving force (Fig. 6). As an additional caveat, an estimation of thermal fatigue in particular rocks with local shadowing may not be a straight-forward extrap-olation of our results and would require using the numerical model with the appropriate temperature inputs.

While our new thermomechanical model constitutes a significant advance from previous semi-analytical models, it is still important to realize the simplifications and inherent assumptions that we have made. For instance, we are now using a two-dimensional mechanical model with thermal inputs from a one-dimensional heat diffusion model. We do not consider the cases where a crack deviates from a straight path, which may cause scabbing of the larger rocks by gradually chipping off material from the surface. A break of the symmetry could occur when using a 2D temperature field, or with the interaction of multiple initial cracks. In these cases, fragments could be produced by crack coalescence or by flaking surface material. Given that an idealized case of a single thermal fatigue crack growth model has demonstrated the importance of this mechanism, it is worthwhile to make future efforts to include the more complex behavior of a network of cracks, as well as the coupling between thermal fatigue and mechanical impact. It is unlikely that rocks on airless bodies would have only a single active surface crack that drives the thermal fragmentation.

Natural geological materials have a large distribution of flaws with different sizes and orientations. These flaws can evolve by thermal fatigue, but are also modified through non-catastrophic micrometeoritic impacts. Impacts by micrometeorites nucleate and propagate cracks, which can then be grown through fatigue. In that sense, these two mechanisms are collaborating processes in asteroidal regolith generation: micro-meteoritic impacts nucleate cracks in rocks that will propagate due to thermal fatigue, which in turn weakens the body making a subsequent impact more effective.

The numerical model developed in this work has the capability to calculate the crack-tip driving forces resulting from the interaction of a large network of cracks. However, it is unclear if the Paris law would still be an appropriate representation of the crack growth rates. The Paris law parameters are empirical fits to fatigue experiments on single crack growth.

In addition, uncertainty quantification through surrogate modeling of fatigue growth [START_REF] Sankararaman | Uncertainty quantification and model validation of fatigue crack growth prediction[END_REF] showed a major contribution of the Paris' C parameter to the errors in predicted fatigue cycles. These parameter uncertainties contribute to large error bars in the fatigue life prediction. [START_REF] Atkinson | Subcritical crack growth in geological materials[END_REF] demonstrated that different forms of the Paris law could yield essentially identical predicted growth rates for laboratory-scale experiments. In geophysical applications, however, we need to extrapolate beyond the bounds of laboratory scales, and different forms of the Paris law diverge substantially, with the form chosen in this work yielding the most conservative estimates. In our simulations, the crack always grew in fatigue mode.

There may be a limit for how small the excursion in stress intensity factors should be for any growth to occur (typically known as the stress intensity threshold). So far, there has not been any clear experimental evidence of the existence of a subcritical crack growth limit or threshold for crack growth in geological material. Indeed, fatigue experiments performed by [START_REF] Wilkins | Slow crack growth and delayed failure of granite[END_REF] on granite showed crack speeds as low as 10 -11 -10 -12 m/s without encountering a crack arrest limit. The reliable prediction of thermal fatigue lifetime of rocks is therefore limited by the absence of a significant database fatigue experiments on planetary materials and in environments comparable to airless bodies (for example, thermal cycling experiments conducted in a vacuum chamber). Regardless of the exact growth law being used (assuming cracks are growing sub-critically in fatigue), we do not expect that the sizedependence we identify in this work to be affected. The crack tip driving forces (shown in Fig. 6) are insensitive to the precise expression of the crack growth law. The rate of breakdown as a function of rock diameter, however, would be greatly affected.

We note that several stress relaxation mechanisms such as creep could occur during these long loading cycles that may modify the stress response and crack growth behavior. These relaxation effects have been reported in sedimentary and hard rocks [START_REF] Maranini | Creep behaviour of a weak rock: experimental characterization[END_REF][START_REF] Fujii | Circumferential strain behavior during creep tests of brittle rocks[END_REF][START_REF] Li | Time-dependent tests on intact rocks in uniaxial compression[END_REF], and can relax peak stresses by 10-15% in roughly 13 hours.

It is not inconceivable, then, that microstructure evolution through creep would relax the induced diurnal thermal stresses over the course of several thousands of years, leading to a slower thermal fatigue cracking.

Mismatches in expansion coefficients between the individual grains could play an important role in driving the thermal fatigue fragmentation of regolith with sizes much smaller than the skin depth. In this work, we consider an averaged contribution to the stress tensor from inclusions by means of a sub-scale representative volume element approach. With this homogenization approach, it is possible that the stresses arising from thermal mismatch alone could be underrepresented, especially in the smaller rocks.

The simplified analytical integration of the Paris law presented in Section 3 showed a good fit to the numerical model results despite using a linear approximation of the excursion stress intensity factors. The derived analyti-cal expressions for the breakdown time can therefore be used for a first-order estimate of thermal fatigue fragmentation for rocks that are of the same order as the skin depth. Rocks with diameters that were much larger than the skin depth showed weak response to thermal cycling, suggesting that the thermal fatigue response would be very weak beyond such sizes. We note also that thermal fatigue scaling is highly nonlinear with changes in rotation period.

This is in a small part due to the longer times associated with each fatigue cycle, but is more dependent on the resulting surface temperature variations and its spatiotemporal evolution. Fast rotators (periods of 2-4 hours) experience more thermal cycles for the same amount of time as slower rotators (periods of 12 hours or more). However, a fast rotator experiences smaller temperature variations per cycle, and therefore smaller driving forces on the crack. This suggests that thermal fatigue would be most efficient on "intermediate" rotators (periods of rotation ∼ 6 hours), and the scaling arguments developed in this work would be best applied to such asteroids. In addition, a layered subsurface in asteroids would greatly affect the temperature profile (and gradients) from that shown in Fig. 1. [START_REF] Mellon | The presence and stability of ground ice in the southern hemisphere of Mars[END_REF] showed that the magnitude of temperature oscillation can be reduced in the subsurface by a factor of 5 in a layered subsurface. Similarly, [START_REF] Molaro | Thermally induced stresses in boulders on airless body surfaces, and implications for rock breakdown[END_REF] showed that the regolith provides a strong insulation to buried rocks on the lunar surface and reduces their temperature excursion by a factor of almost 18. For these cases, which in our study would be analogous to a reduction in both the skin depth and the temperature excursion, we would not expect to see significant thermal fatigue crack growth in those buried rocks.

Modeling thermal fatigue is a computationally intensive process that re-quires careful tracking of the crack front over a large number of cycles. The computational framework needs to be able to bridge efficiently between the relatively fast rotational periods (several hours), to the final fragmentation time of a rock (thousands to millions of years). The numerical model presented here is designed to efficiently model the fatigue growth of cracks in regolith over a several thousands of years, while using timesteps as small as 15 minutes. More sophisticated fracture codes, such as the mapped finite element method [START_REF] Chiaramonte | Mapped finite element methods: High-order approximations of problems on domains with cracks and corners[END_REF], have been developed recently.

The mapped finite element method provides a high-order approximation of problems with cracks and yields optimal convergence rates without additional degrees of freedom or special shape functions. Such models may provide improved accuracy at reduced computational costs and their implementation in the code developed in this work could be potential paths for improving the numerical algorithm. This advance in numerical modeling capabilities further emphasizes the need for careful baseline thermal fatigue experiments that constrain the material parameters and crack growth laws that are used in such predictive models.

Conclusions

We find that the rate of rock breakdown by thermal fatigue has two domains that are described by distinct slopes (Fig. 5). Using an approximate integration of the Paris law based on the results of the detailed fracture simulations, we derive measures for these two slopes. Normalized plots of computed fragmentation times for different rock diameters seem to align well under these analytically derived curves. Our analysis demonstrates that

Figure 1 :

 1 Figure 1: Temperature in a rock at 1 AU as a function of depth during a 6 hour day/night cycle, plotted at 15 min increments. This temperature profile is typical of diurnal thermal cycling. The heat wave attenuates quickly beyond the skin depth, and temperature variations become negligible.

Figure 3 :

 3 Figure 3: Flowchart of the thermomechanical model developed in this work. The preprocessing stage sets up the rock geometry, thermomechanical properties, and initial conditions. The processing stage runs in parallel on multiple computational nodes to solve different timesteps within a single period of rotation. The excursion in stress intensity factor is calculated in the post-processing stage and used to advance cracks until fragmentation.

  presents the computed stress field in a 40 cm diameter rock containing a 15 cm crack. The figure shows six snapshots (the columns) from the computations, one every hour over the 6 hour period. The top row shows the temperature field in the rock at that time, and the bottom row shows the corresponding stress field in the cracked rock. For comparison, the middle row shows the stress field that would exist in the uncracked rock. The cracked rock has an entirely different stress field than the uncracked rock, with different maximum stresses, and with these maxima located in different parts of the rock. These differences highlight the importance of using fracture-mechanics solutions to accurately capture the singularity at the crack tip. Note that answering the question of how much a crack will grow in such a thermal cycle can only be done by using the stress intensity factor as a measure. Further, stresses alone cannot identify the extent of crack growth. Thus the peak stresses obtained from an uncracked rock alone are not reasonable indicators of the thermal fatigue lifespan of cracked rocks.

Figure 4 :

 4 Figure4: The temperature field in a 40 cm diameter rock is plotted in the first row at 1 hour increments for a 6 hour thermal cycle. The stress field perpendicular to the crack faces (stresses in the crack opening direction) is shown on the second and third rows. The second row is for a rock without any crack, while the third row shows the stress field in a rock of the same size that has a 15 cm surface crack. Notice that the stress at the crack tip is singular and that the crack relaxes the stress field at its faces.

  thermal fatigue experiments that have been conducted on meteorite samples are those conducted by[START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF] and[START_REF] Hazeli | The origins of asteroidal rock disaggregation: Interplay of thermal fatigue and microstructure[END_REF]. The observed crack growth in the experiments of[START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF] is well represented by a thermo-mechanical model similar to the one described here and a crack growth rate following the Paris law with parameters values similar to those of Carrara marble (from[START_REF] Migliazza | Experimental investigation on crack propagation in Carrara marble subjected to cyclic loads[END_REF]).[START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF]

  Fig. 5. The figure shows the time taken (in years) to fragment a rock of a given size (the blue data points are obtained from our simulations). All rocks are assumed to have an initial 30µ m crack. Due to the numerical resolution and the nature of the contour integral, this crack typically grows to around 0.1 of the rock diameter during the first few cycles. As such, we do not expect a large deviation in the survival times as a result of the initial crack size. The bars around each data point represent a variation of 30% in the Paris coefficientC. Note that the computed survival times can be as high as 10 8 years, while the thermal cycling period is 6 hours, demonstrating the need for an efficient computational scheme to handle this wide range of timescales. For rocks smaller than 7 cm, our model predicts rock survival times (Fig.5, in blue) that are consistent with those reported in[START_REF] Delbo | Thermal fatigue as the origin of regolith on small asteroids[END_REF] (Fig.5, in red) using lower-accuracy computations. We predict longer fragmentation times thanDelbo et al. (

Figure 5 :

 5 Figure 5: Survival time of cm-sized (diameter) carbonaceous chondrites at 1 AU and a dirunal cycle of 6 hours. The data points in red are from Delbo et al. (2014), while the blue data points are obtained from the model presented in this work. Error bars represent uncertainties (30% variation) in the Paris' parameters. The dashed line shows the times at which 50% of these same rocks are broken by micrometeoroid impacts, using the Hörz et al. (1975) model with the appropriate modifications for asteroids. The survival times of rocks smaller than the diurnal skin depth (approx. 7 cm) follows the same trend as reported in Delbo et al. (2014): larger rocks are fragmented faster than smaller rocks. However, as rock diameters become larger than the skin depth, a reversal in the trend is observed. That is, the crack growth speed of thermal fatigue driven surface crack is greatly reduced in these larger rocks.

Figure 6 :

 6 Figure6: Normalized excursion in stress intensity factor as a function of crack length normalized by the rock diameter. In small (less than 10 cm) rocks, ∆K/K cr generally increases with increasing crack size, especially between a/d of 0.1 and 0.7. Rocks much larger than the skin depth see a different trend, wherein the excursion in stress intensity factor decreases as the crack size increases. The implications of these trends on the survival rates of asteroid rocks are discussed in Section 3.

Figure 7 :

 7 Figure7: Scaled survival times of carbonaceous chondrites (CC) and ordinary chondrites (OC) on near-Earth asteroids and in the main belt. Circles and squares represent numerical simulation results for main belt and Near-Earth asteroids, respectively. The x-axis is a measure of the survival rates normalized by the number of cycles for fragmenting a rock with diameter equal to the skin depth. It is defined from Eq. (23) as t = (t f )/t δf , where t δf is the fragmentation time for a rock with diameter equal to the skin depth. The y-axis is the rock diameter normalized by the skin depth. With this dimensionless representation, the two slopes intersect at (1, 1). The circle and square symbols are results from the numerical model. The solid lines represent the slopes identified in Eq. (16a) and the shaded region represents 30% confidence intervals related to the Paris' exponent n.

Figure 8 :

 8 Figure8: Predicted time (in years) to break down a 10 cm diameter rock on a C-type asteroid as a function of period of rotation and heliocentric distance, using the scaling factors in Eq. (24). Some C-type asteroids are highlighted on the plot using their semi-major axis and rotation period as reference. These were not modeled directly, and discrepancies from Table1in their particular surface and material composition would lead to different breakdown rates than those predicted in this figure. These results indicate the existence of an "intermediate" range in rotational periods where thermal fatigue would be more efficient. Very fast rotators have small skin depths, and larger rocks would require increasingly larger times to fragment. Slow rotators have larger skin depths and take longer time to complete a single period, and rocks much smaller than that skin depth would also take increasingly longer times to fragment (see discussion for Fig.5). Consequently, the fragmentation times and the ideal "intermediate range" in rotation period depend on the particular rock diameter and how it compares with the skin depth.

Table 1 :

 1 Values and base dimensions assumed for the physical properties of the carbonaceous chondrite material considered in this study. The base dimensions are expressed in the MLTθ (Mass-Length-Time-Temperature) system.

		Symbol Units			Dimensions	Carbonaceous Chondrite Ordinary Chondrite	Reference
	Paris exponent Paris pre-factor	n C	-m (MPa	√	m) -n	-L 1+n/2 M -n T 2n 3 × 10 -4 3.84	3.84 3 × 10 -4	Migliazza et al. (2011) Migliazza et al. (2011)
	Young's Modulus	E	GPa			M L -1 T -2	45	74	Delbo et al. (2014)
	Poisson's ratio	ν	-			-	0.24	0.28	Delbo et al. (2014)
	Bulk expansion coefficient	αm	K -1			θ -1	8.5 × 10 -6	8.5 × 10 -6	Hazen (1977)
	Chondrule expansion coefficient αinc	K -1			θ -1	10.4 × 10 -6	10.4 × 10 -6	Smyth (1975)
	Thermal Conductivity	κ

rocks with diameters much larger than the diurnal skin depth will exhibit poor thermal fatigue cracking. In such rocks, cracks can grow quickly by thermal fatigue up to a certain length, after which their progress is slowed down greatly (or even arrested completely).

The identification of this critical lengthscale provides bounds on the areas where thermal fatigue is likely to be the dominant mechanism in rock disaggregation. We expect that large boulders would not be fragmented completely by thermal fatigue, but would have long straight cracks that grow in the direction of maximum circumferential thermal stress, just as in the Martian boulders shown in [START_REF] Eppes | Cracks in Martian boulders exhibit preferred orientations that point to solarinduced thermal stress[END_REF]. We expect that this critical lengthscale may be an important indicator for the expected regolith size distribution. The thermal fatigue contribution may thus be reflected as deviations in the regolith particle size-frequency distribution as a result of preferentially depleting rocks that are of dimensions comparable to the skin depth.

Acknowledgements