Mikael Beuthe
email: mbeuthe@oma.be

East-west faults due to planetary contraction

Keywords:

Contraction, expansion and despinning have been common in the past evolution of solar system bodies. These processes deform the lithosphere until it breaks along faults. Their characteristic tectonic patterns have thus been sought for on all planets and large satellites with an ancient surface. While the

search for despinning tectonics has not been conclusive, there is good observational evidence on several bodies for the global faulting pattern associated with contraction or expansion, though the pattern is seldom isotropic as predicted. The cause of the non-random orientation of the faults has been attributed either to regional stresses or to the combined action of contraction/expansion with another deformation (despinning, tidal deformation, reorientation). Another cause of the mismatch may be the neglect of the lithospheric thinning at the equator or at the poles due either to latitudinal variation in solar insolation or to localized tidal dissipation. Using thin elastic shells with variable thickness, I show that the equatorial thinning of the lithosphere transforms the homogeneous and isotropic fault pattern caused by contraction/expansion into a pattern of faults striking east-west, preferably formed in the equatorial region. By contrast, lithospheric thickness variations only weakly affect the despinning faulting pattern consisting of equatorial strike-slip faults and polar normal faults. If contraction is added to despinning, the despinning pattern first shifts to thrust faults striking north-south and then to thrust faults striking east-west. If the lithosphere is thinner at the poles, the tectonic pattern caused by contraction/expansion consists of faults striking north/south. I start by predicting the main characteristics of the stress pattern with symmetry arguments. I further prove that the solutions for contraction and despinning are dual if the inverse elastic thickness is limited to harmonic degree two, making it easy to determine fault orientation for combined contraction and despinning. I give two methods for solving the equations of elasticity, one numerical and the other semi-analytical.

The latter method yields explicit formulas for stresses as expansions in Legendre polynomials about the solution for constant shell thickness. Though I only discuss the cases of a lithosphere thinner at the equator or at the poles, the method is applicable for any latitudinal variation of the lithospheric thickness. On Iapetus, contraction or expansion on a lithosphere thinner at the equator explains the location and orientation of the equatorial ridge. On Mercury, the combination of contraction and despinning makes possible the existence of zonal provinces of thrust faults differing in orientation (north-south or east-west), which may be relevant to the orientation of lobate scarps.

Introduction

A great variety of tectonic features is found on nearly all solid planets and large satellites of the solar system: ridges and scarps, rifts and grabens, furrows and grooves etc. Their origin on Earth mainly lies in the movement of nearly rigid plates but other mechanisms must be found elsewhere since plate tectonics is unique to Earth. Though tectonic features can often be explained by a regional effect, such as an impact or the emplacement of a volcanic load, a subclass of them sometimes form a global pattern on the surface. In such a case, their cause must be the global deformation of the whole lithosphere which generates a global stress distribution resulting in a characteristic faulting pattern at the surface.

Several processes affect the planetary figure. Though a change in the mean planetary radius (contraction or expansion) is the simplest conceivable deformation, its underlying causes can be complex since they depend on the complicated physics of the interior of a self-gravitating differentiated body (e.g. [START_REF] Andrews-Hanna | Strike-slip faults on Mars: Observations and implications for global tectonics and geodynamics[END_REF] for Mars, [START_REF] Squyres | The tectonics of icy satellites[END_REF] and [START_REF] Collins | Tectonics of the Outer Planet Satellites[END_REF] for icy satellites). On a planet with a lithosphere of constant thickness, contraction and expansion bring about a homogeneous distribution of randomly oriented compressional and extensional faults, respectively [START_REF] Melosh | Global tectonics of a despun planet[END_REF]. The surfaces of Mars [START_REF] Anderson | Primary centers and secondary concentrations of tectonic activity through time in the western hemisphere of Mars[END_REF][START_REF] Anderson | Centers of tectonic activity in the eastern hemisphere of Mars[END_REF][START_REF] Knapmeyer | Working models for spatial distribution and level of Mars' seismicity[END_REF] and Mercury (Watters et al., 2009) show widespread compressional tectonic features, termed wrinkle ridges and lobate scarps. Their distribution is however far from uniform and their orientation is not random. On Mars, regional effects such as the Tharsis rise have strongly influenced the global pattern of wrinkle ridges and the possible role of a global contraction event remains under discussion [START_REF] Mangold | Chronology of compressional deformation on Mars: evidence for a single and global origin[END_REF][START_REF] Nahm | A test of global contraction models for Mars using observations[END_REF]. On Mercury, the lobate scarps that have been identified in Mariner 10 images have a greater cumulative length in the southern latitudes, generally trend within 50 • of the north-south direction and preferably dip northward below 50 • S [START_REF] Watters | Thrust faults and the global contraction of Mercury[END_REF][START_REF] Watters | The Tectonics of Mercury[END_REF]. The anisotropy of lobate scarps is often attributed to the additional contribution of despinning (see below), but the regional effect of the Caloris basin has also been invoked [START_REF] Thomas | Tectonic history of Mercury[END_REF]. A more complete and less illumination-biased fault catalog should result from the analysis of Messenger images and altimetry data; results from the first flyby are promising [START_REF] Solomon | Return to Mercury: A Global Perspective on MESSENGER's First Mercury Flyby[END_REF]Watters et al., 2009). Among large icy satellites, Ganymede is a showcase for extensional tectonics [START_REF] Pappalardo | Geology of Ganymede[END_REF].

Global expansion is also thought to have played a role in the formation of the global tectonic grids on Rhea [START_REF] Consolmagno | Resurfacing Saturn's satellites -Models of partial differentiation and expansion[END_REF][START_REF] Moore | The geomorphology of Rhea: implications for geologic history and surface processes[END_REF][START_REF] Thomas | The tectonic and volcanic history of Rhea as inferred from studies of scarps, ridges, troughs, and other lineaments[END_REF], Dione [START_REF] Moore | The tectonic and volcanic history of Dione[END_REF][START_REF] Consolmagno | Resurfacing Saturn's satellites -Models of partial differentiation and expansion[END_REF] and Ariel [START_REF] Plescia | Geological terrains and crater frequencies on Ariel[END_REF][START_REF] Nyffenegger | Tectonic lineations and frictional faulting on a relatively simple body (Ariel)[END_REF], though the orientation of faults has been shown to be non-random in each case. Global contraction may also have occurred on Rhea and Dione. New tectonic analyses of Rhea [START_REF] Moore | Topography of Endogenic Features on Saturnian Mid-Sized Satellites[END_REF][START_REF] Wagner | The Global Geology of Rhea: Preliminary Implications from the Cassini ISS Data[END_REF] and Dione [START_REF] Goff-Pochat | Strain Measurement Across Fault Scarps on Dione[END_REF][START_REF] Wagner | Stratigraphy of Tectonic Features on Saturn's Satellite Dione Derived from Cassini ISS Camera Data[END_REF][START_REF] Stephan | Dione's spectral and geological properties[END_REF] are underway with Cassini data. While partial imaging makes it difficult to identify a global tectonic grid on Iapetus [START_REF] Singer | A Search for Despinning Fractures on Iapetus[END_REF], its huge equatorial ridge must be related to a global deformation.

Contraction [START_REF] Castillo-Rogez | Iapetus' geophysics: Rotation rate, shape, and equatorial ridge[END_REF] has been cited as a possible culprit, but the corresponding faulting pattern neither predicts the equatorial location nor the east-west orientation of the ridge.

The next simplest deformation is a change in the planetary flattening due to a decrease in the rotation rate, or despinning, that results from tidal effects leading to a resonant or synchronous rotation [START_REF] Murray | Solar System Dynamics[END_REF]. Despinning tectonics on a planet with a thin lithosphere of constant thickness consist of an equatorial zone of strike-slip faults (striking at about 60 degrees from the north) and of polar zones of extensional faults or joints striking east-west [START_REF] Burns | Consequences of the tidal slowing of Mercury[END_REF][START_REF] Melosh | Global tectonics of a despun planet[END_REF]. Despun bodies are common in the solar system:

Mercury is in resonant rotation with the Sun [START_REF] Goldreich | The Dynamics of Planetary Rotations[END_REF]) and nearly all large satellites in the solar system are in synchronous rotation with their parent body [START_REF] Peale | Origin and Evolution of the Natural Satellites[END_REF]. Simultaneous despinning and contraction have been used to justify the predominant north-south orientation of lobate scarps on Mercury [START_REF] Melosh | Mercurian global tectonics -A consequence of tidal despinning[END_REF][START_REF] Pechmann | Global fracture patterns of a despun planet -Application to Mercury[END_REF][START_REF] Dombard | Despinning plus global contraction and the orientation of lobate scarps on Mercury: Predictions for MESSENGER[END_REF]. Beside the pattern of young lobate scarps, Mercury exhibits a global grid of more ancient lineaments which [START_REF] Melosh | Mercurian global tectonics -A consequence of tidal despinning[END_REF] and [START_REF] Pechmann | Global fracture patterns of a despun planet -Application to Mercury[END_REF] interpreted as evidence of despinning (plus contraction), though the case is not closed [START_REF] Melosh | The tectonics of Mercury[END_REF][START_REF] Thomas | Tectonic history of Mercury[END_REF]. While despinning is a generic phenomenon for large icy satellites, none exhibits unambiguous evidence of the corresponding global tectonic pattern, either because of later resurfacing or because despinning occurred before the surface was fully formed [START_REF] Squyres | The tectonics of icy satellites[END_REF]. The location of Iapetus' ridge suggests a relation with despinning [START_REF] Porco | Cassini Imaging Science: Initial Results on Phoebe and Iapetus[END_REF][START_REF] Castillo-Rogez | Iapetus' geophysics: Rotation rate, shape, and equatorial ridge[END_REF] but it was immediately noted that despinning tectonics cannot produce east-west features at the equator. Global tectonic patterns can also be produced by other types of deformations not examined in this paper: reori-entation relative to the spin axis [START_REF] Melosh | Tectonic patterns on a reoriented planet -Mars[END_REF][START_REF] Leith | Is There Evidence for Polar Wander on Europa?[END_REF][START_REF] Matsuyama | Tectonic patterns on reoriented and despun planetary bodies[END_REF], recession (Melosh, 1980a;[START_REF] Helfenstein | Patterns of fracture and tidal stresses on Europa[END_REF], diurnal tides and non-synchronous rotation [START_REF] Helfenstein | Patterns of fracture and tidal stresses due to nonsynchronous rotation -Implications for fracturing on Europa[END_REF][START_REF] Greenberg | Tectonic Processes on Europa: Tidal Stresses, Mechanical Response, and Visible Features[END_REF][START_REF] Wahr | Modeling stresses on satellites due to nonsynchronous rotation and orbital eccentricity using gravitational potential theory[END_REF]. These deformations can be seen as a superposition of biaxial deformations and thus share the basics of the despinning model.

Up to now the distribution and orientation of tectonic features have always been computed with the assumption of a lithosphere of constant thickness, one reason being the lack of data about the lithospheric thickness and its variations. The choice of a specific variation of lithospheric thickness may thus seem ad hoc if its cause is a variation in the internal heating flux. Two external phenomena provide a better motivation for a model of lithospheric thickness variation. On the one hand, the latitudinal variation in the solar insolation elevates the lithospheric isotherm at the equator of a planet without atmosphere.

For Mercury, the lithosphere could be thinner by a factor of two at the equator in comparison with the poles [START_REF] Mckinnon | Application of ring tectonic theory to Mercury and other solar system bodies[END_REF][START_REF] Melosh | The tectonics of Mercury[END_REF]; a longitude variation is also expected because of the difference in average temperature between the so-called equatorial hot and cold poles. On the other hand, tidal heating acting on the satellites with an eccentric orbit creates additional lateral variations in the lithospheric thickness. [START_REF] Nimmo | The global shape of Europa: Constraints on lateral shell thickness variations[END_REF] however found no evidence in the global shape of Europa for the thickness variations predicted by [START_REF] Ojakangas | Thermal state of an ice shell on Europa[END_REF] (see also [START_REF] Tobie | Tidally heated convection: Constraints on Europa's ice shell thickness[END_REF]). Since large satellites are at present in a synchronous state of rotation, tidal heating is localized in longitude. Nevertheless tidal heating is evenly distributed in longitude before tidal locking is achieved, so that the lithospheric thickness only depends on latitude during the despinning phase.

The goal of this paper is to predict global tectonic patterns due to contraction and despinning on a lithosphere with latitudinal thickness variation. Though I will assume in the numerical examples that the thickness variation is symmetrical about the equatorial plane, the method is also valid for non-symmetrical variations. The prediction of tectonic patterns follows the following procedure.

First the lithosphere is modeled as a thin elastic spherical shell [START_REF] Turcotte | Role of membrane stresses in the support of planetary topography[END_REF][START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]. Second the contraction or despinning is represented as a static deformation of the shell given by the difference between the initial and final figure of the planet. The stresses caused by the deformation of the thin elastic spherical shell are computed with the equations of equilibrium and elasticity. Third the style and orientation of faults is predicted from the stresses with Anderson's theory of faulting (Jaeger et al., 2007).

The most significant result of this paper concerns the contraction of a planet with a lithosphere thinner at the equator (predictions for expansion follow by changing the sign of the stresses). If the contracting lithosphere is of constant thickness, stresses are homogeneous and isotropic, which means that their orientation and magnitude are independent of position and that the principal tangential stresses are equal. The thinning of the lithosphere at the equator has two effects: (1) the stress becomes most compressive at the equator and (2) the meridional component of the stress becomes more compressive than the azimuthal component. This situation, as well as the corresponding expansion case, favors the development of faults striking east-west and starting preferably at the equator. The contraction of a shell thinner at the equator is thus complementary to despinning, for which the azimuthal stress is always more compressive than the meridional stress. If the inverse thickness is at most of harmonic degree two I show that this complementarity is described by an exact relation between a contracting planet and a planet that is spinning up. As for despinning, its tectonic pattern is not much affected by the variation of the lithospheric thickness: the main effect is the reduction in size of the polar provinces of tensional stress as the lithosphere thickens at the poles. The combination of contraction and despinning renders possible the existence of provinces of thrust faults differing in their orientation (north-south or east-west), but this only occurs for a specific ratio between contraction and despinning. If the lithosphere is thinner at the poles, global contraction (resp. expansion) leads to thrust (resp. normal) faults striking north-south and preferably formed in the polar regions, whereas the despinning pattern is weakly affected.

As for applications, I show that contraction or expansion combined with lithospheric thinning at the equator provides an explanation for the location and orientation of Iapetus' equatorial ridge. I also attempt to account for the orientation of lobate scarps on Mercury. Tectonic patterns resulting from simultaneous contraction and despinning on a lithosphere thinner at the equator are consistent with the latitudinal dependence of lobate scarps. The amount of required contraction is however much larger than what has been estimated from observations, making it necessary to resort to more contrived scenarios involving fault reactivation.

The paper can be read at two levels: I present general results and geophysical interpretations in the main text, while the semi-analytical method of solution is described in detail in the appendices, which also include explicit solutions as well as a discussion of the relationship between thin shells, thick shells and Love numbers. Section 2 states the membrane equations that must be solved in order to compute the stresses generated by the deformation of a thin elastic shell with variable thickness. Remarkably, assumptions of axial and equatorial symmetry (presented in Section 3) are sufficient for the prediction of the basic characteristics of stress and strain. These symmetries lead to dualities between the contraction and despinning solutions when the inverse thickness is at most of degree two, reducing by half the computational load and simplifying the discussion of combined contraction and despinning. I also compute a first order approximation of the contraction stresses from duality arguments. In Section 4, I show how to solve the membrane equation. After discussing the parameterization of the thickness variation, I frame the problem for a numerical solution with Mathematica, before explaining a semi-analytical method in which truncated Legendre expansions serve to rewrite the membrane equation as a matrix equation. In Section 5, I use the contraction and despinning solutions of the membrane equation in order to predict the extent of tectonic provinces, i.e. areas with faults of given style and orientation. The cases of expansion and spinning-up are related to contraction and despinning by a sign change. In Section 6, I discuss the application of the predicted tectonic patterns to Iapetus' ridge and Mercury's lobate scarps.

Membrane stress in a deformed spherical shell

The style and orientation of faults due to contraction and despinning can be predicted from knowledge of the stresses in the lithosphere. The lithosphere is modeled as a thin elastic spherical shell (in this paper, lithospheric thickness and elastic thickness are synonymous). The important assumptions that the shell is thin and that it is elastic are analyzed by [START_REF] Melosh | Global tectonics of a despun planet[END_REF] for a shell of constant thickness.

First, Melosh shows that models with a thin or a thick shell predict very similar despinning faulting patterns, except the possible presence of an equatorial thrust fault province when the shell is thick (the contraction tectonic pattern is not affected by the thickness of the shell). Despinning stresses at the surface of a thick shell can easily be computed if the secular tidal Love numbers h T 2 and l T 2 are known (see D). At a given latitude, the stress magnitude mainly depends on h T 2 whereas the relative size of the principal stresses depends on the ratio l T 2 /h T 2 . When the shell thickness increases, h T 2 decreases so that despinning stresses decrease in magnitude. While Melosh computes thin shell stresses under the assumption of hydrostatic flattening, [START_REF] Matsuyama | Tectonic patterns on reoriented and despun planetary bodies[END_REF] extend the domain of application of thin shell formulas by using a nonhydrostatic flattening parameterized by h T 2 . This latter procedure is equivalent to keeping the overall h T 2 -dependence in thick shell formulas while setting the ratio l T 2 /h T 2 to its membrane value (1 + ν)/(5 + ν), where ν is Poisson's ratio (see D).

Second, Melosh argues that the faulting pattern predicted by the elastic model is not substantially modified when the lithosphere becomes plastic, though the numerical values of the stresses are bounded by the yield stress. The elastic model is thus limited to the prediction of the faulting style (including its orientation) but does not say whether faulting occurs.

Stresses in a shell satisfy equilibrium equations which relate them to the load deforming the shell. However the load causing the contraction or the flattening change during despinning is a priori unknown. Thus the full formulation of elastic equations involving stress-strain and strain-displacement relations becomes necessary. For an elastic shell of constant thickness, contraction stresses were obtained by Lamé (see Love, 1944, p. 142) whereas despinning stresses were computed by [START_REF] Vening-Meinesz | Shear patterns of the Earth's crust[END_REF]. There are no analytic formulas giving equivalent results for a thin elastic shell with variable thickness but [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF] recently derived the scalar equations governing the deformation of such a shell under an arbitrary load.

Since displacements of the shell surface due to contraction and despinning have a large wavelength, bending moments in the shell are negligible: the shell mainly responds by stretching and is said to be in the membrane regime [START_REF] Kraus | Thin Elastic Shells[END_REF][START_REF] Turcotte | Role of membrane stresses in the support of planetary topography[END_REF]. In the membrane limit, bending moments are set to zero: deformation equations take a simplified form which can be obtained from the full deformation equations by setting to zero the bending rigidity D. In that case, the membrane equations for a spherical shell of radius R and variable thickness h are ∆ F = -R q ,

(1)

C (α ; F ) -(1 + ν) A (α ; F ) = 1 R ∆ w . ( 2
)
where q is the transverse load (positive inward), w is the transverse displacement (positive outward) and F is an auxiliary function called the stress function.

Tensile stress is positive. The tangential load has been set to zero. The differential operators ∆ , C and A are defined in A by Eqs. ( 63)-( 67). The thickness h is included in the elastic parameter α:

α = 1 Eh . ( 3
)
E is Young's modulus (which is allowed to be variable), ν is Poisson's ratio (which must be constant).

In the membrane limit, the stress integrated over the thickness of the shell, or stress resultant, is given by

(N θθ , N ϕϕ , N θϕ ) = (O 2 , O 1 , -O 3 ) F . ( 4
)
with the operators O i being defined by Eqs. ( 60)-( 62). The stress as expressed in Eq. ( 4) automatically satisfies the two tangential equilibrium equations while the radial equilibrium equation is encoded in Eq. (1). In the membrane limit, there is thus only one equation where elastic parameters appear, Eq. ( 2), which results from an equation of compatibility between the strains.

In the membrane limit, the stress is constant through the shell thickness. It is thus related to the stress resultant by

σ ij = 1 h N ij . ( 5
)
The strain is given by

θθ = 1 E (σ θθ -νσ ϕϕ ) , ( 6
) ϕϕ = 1 E (σ ϕϕ -νσ θθ ) , ( 7
) θϕ = 1 + ν E σ θϕ . ( 8
)
If the load is known, the stress function F can be computed from the first membrane equation (1) and the stress from Eqs. (4)-( 5). In that case, the elastic thickness only appears in Eq. ( 5), which means that the solution for variable thickness is directly obtained from the solution of constant thickness by the scaling h 0 /h (h 0 is the mean elastic thickness). This procedure implies that the load is the same whether the elastic thickness is constant or not. However the load is a priori unknown when modeling contracting and despinning events; instead, the transverse displacement w is the known input. The stress function F must thus be computed from the second membrane equation, Eq. ( 2), which is linear in F and w with variable coefficients depending on 1/Eh. As a corollary, the load causing contraction or despinning is not the same in the cases of constant and variable elastic thickness.

It should be noted that the variation of the lithospheric thickness has an effect on the deformation. Contraction is not a problem if it is defined as a uniform change in radius, but it is another matter if the definition involves physical processes such as an internal thermal contraction or an expansion due to water-ice transition. Despinning causes deformations of harmonic degrees zero and two (the former is very small and usually neglected) if the lithospheric thickness is constant but other harmonic degrees appear if it is variable. In principle the effect of internal contraction and despinning on the shape of a planet with a variable lithospheric thickness should be computed with a threedimensional model of the interior. Love numbers quantify the response of a planet with a spherically symmetric internal structure (see D), but there is no standard technique dealing with the non-spherical case. Extreme cases of a very thin or a very thick lithosphere are not a problem. On the one hand, the lithosphere does not affect the shape of the body in the limit of vanishing thickness. On the other hand, thickness variations of a thick lithosphere are expected to be relatively small with respect to the average thickness: latitudinal variations in solar insolation, for example, induce smaller thickness variations (in percentage) in a thick lithosphere than in a thin one. In both cases, the response of the planet is thus adequately described by using a model with a spherically symmetric internal structure. Between these extremes, the size of the deformations with harmonic degrees not equal to zero (for contraction) or two (for despinning) should be investigated, but this question is beyond the scope of this paper. I thus choose to define contraction and despinning as phenomena associated with deformations of harmonic degree zero and two, respectively.

How good an approximation is the membrane limit? If bending terms are not neglected, Eq. ( 2) remains the same but new terms appear in Eq. ( 4) for the stress resultants. These new terms have the form (D/R 3 ) O i w (with D = Eh 3 /12(1ν 2 )) so that they are smaller than the dominant term O i F by a factor ξ = 12R 2 /h 2 . Neglected bending terms are thus of the same order of magnitude as other small terms present in the equations before the 'thin shell' limit of large ξ is taken [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]. It is possible to solve the complete equations (with bending terms and no large ξ limit) with the same method as proposed in this paper, but the gain in precision would only be of order 1/ξ.

Even so the theory remains a thin shell theory in the sense that the radial stress is neglected.

The membrane and stress equations ( 2)-( 5) become nondimensional with the following definitions:

ᾱ = α α 0 , ( 9
) F = α 0 F , ( 10
) w = w R , ( 11
) Nij = α 0 N ij , ( 12
) σij = 1 E σ ij , ( 13
)
where α 0 is the average of α on the sphere, i.e. the coefficient of degree zero in the expansion of α in Legendre polynomials (α 0 = 1 Eh0 ). The membrane equation (2) then becomes

C ᾱ ; F -(1 + ν) A ᾱ ; F = ∆ w . ( 14
)
The nondimensional stress is computed with the nondimensional version of Eqs. (4)-( 5):

Nθθ , Nϕϕ , Nθϕ = (O 2 , O 1 , -O 3 ) F , ( 15
) σij = ᾱ Nij . ( 16
)
The strain is already nondimensional:

θθ = ᾱ (O 2 -νO 1 ) F , ( 17
) ϕϕ = ᾱ (O 1 -νO 2 ) F , ( 18
) θϕ = -ᾱ (1 + ν) O 3 F . ( 19
)
The solution of the nondimensional membrane equation ( 14) when w is only of degree zero ( w = w0 P 0 ) will be called the contraction solution F C (whatever the sign of w0 ), whereas the solution when w is only of degree two ( w = w2 P 2 ) will be called the despinning solution F D (whatever the sign of w2 ). Negative (resp.

positive) values of w0 correspond to contracting (resp. expanding) planets, whereas positive (resp. negative) values of w2 correspond to despinning (resp.

spinning-up) planets. Solutions for a constant elastic thickness are computed in C.

Symmetry and duality

As discussed in the introduction, I make the simplifying assumption that the value of the elastic thickness only depends on the latitude. The contraction and despinning problems are then axially symmetric. Though the numerical method developed in Section 4 does not require it, I also assume equatorial symmetry, i.e. mirror symmetry about the equatorial plane. These two assumptions approximately hold for the physical causes of the thickness variation considered here: variation in solar insolation and tidal heating during despinning, provided that the obliquity of the body with respect to the Sun (in the former case) or to the tidal source (in the latter case) is small.

Using symmetry and duality arguments, I describe in this section the characteristics of stress and strain on contracting or despinning shells thinner (or thicker) at the equator. The duality relation will also serve to obtain a numerical approximation of the contraction solution.

Axial and equatorial symmetry

The main characteristics of the contraction solution can be determined from symmetry arguments, without solving the membrane equation. Axial symmetry give the following constraints on stress and strain (see also Fig. 1):

1. the meridional and azimuthal stresses (resp. strains) are principal stresses (resp. strains), i.e. σ θϕ = 0 (resp. θϕ = 0).

2. the meridional and azimuthal stresses are equal at the poles; this is also true of the strains.

3. the slopes of the stress and strain components vanish at the poles (it is a consequence of the previous point combined with the equations of equilibrium); they also vanish at the equator if there is equatorial symmetry.

4. the average of the meridional strain is independent of the elastic thickness (see Eq. ( 116)): < θθ >= w0 for contraction, < θθ >= w2 /4 for despinning.

5. the azimuthal strain at the equator is independent of the elastic thickness (see Eq. ( 118)): ϕϕ | θ= π 2 = w0 for contraction, ϕϕ | θ= π 2 = -w2 /2 for despinning.

[Figure 1 about here.] These constraints allow me to predict the response of the shell to a radius variation. The stresses for a constant elastic thickness are given in C. Fig. 1 illustrates the following discussion. I suppose that the shell is expanding ( w0 > 0) so that stresses and strains are positive (the solution for a contracting shell is obtained by changing the sign of stresses and strains). I also assume equatorial symmetry with a monotonous variation of the thickness from the pole to the equator, so that stress and strain components can be expected to vary monotonously. Stress and strain are concentrated where the shell is weaker. In other words, the membrane stretching is maximum where the shell is thinnest. If the shell is thinner at the equator, the azimuthal strain is thus maximum at the equator, where it takes the value w0 imposed by the equatorial length variation, and minimum at the poles, where it is equal to the meridional strain. The meridional strain is also maximum at the equator and minimum at the poles, taking somewhere in between its average value w0 . Therefore, the meridional strain always exceeds the azimuthal strain when the shell is thinner at the equator. This ordering is conserved for the stresses: the meridional stress always exceeds the azimuthal stress when the shell is thinner at the equator, as shown by the inversion of Eqs. ( 6)-(8) with ν < 1. If the shell is thinner at the poles, a similar analysis leads to the conclusion that the strain (resp. stress) is maximum at the poles and that the azimuthal strain (resp. stress) always exceeds the meridional strain (resp. stress).

The despinning shell case ( w2 > 0) is not as simple since the combined effects of the non-spherical deformation and thickness variation can lead to local extrema of the stress and strain components between the poles and the equator.

However the ordering of the strains ( ϕϕ more compressive than θθ ) is expected to be the same as in case of constant elastic thickness, since the meridional strain has an average positive value ( w2 /4) while the azimuthal strain varies between the value of the meridional strain at the poles and the negative value -w2 /2 at the equator. The location of the shell thinning has thus less influence on the stress/strain pattern than in the contraction case. Because of stress/strain concentration in weaker areas, the stress and strain components are expected to come closer to (resp. farther from) zero at the poles if the shell becomes thinner at the equator (resp. at the poles), as shown in Fig. 1b. The meridional strain curve will accordingly move so as to keep its average constant. The principle of stress/strain concentration does not apply to the azimuthal strain at the equator since its value is fixed by the geometry. As above, the ordering is conserved when going from strains to stresses: σ ϕϕ is more compressive than σ θθ .

Another argument for ordering will be obtained from the contraction-despinning duality in Section 3.2.

Contraction-despinning duality

Still assuming axial and equatorial symmetry, I now show that the contraction and despinning solutions are not independent if the harmonic expansion of the inverse thickness is limited to degree two. Even if this assumption is not strictly valid, the harmonic of degree two will dominate if the thickness variation is monotonous between the poles and the equator. Under these assumptions, ᾱ can be written in terms of Legendre polynomials of degrees zero and two:

ᾱ = 1 + ᾱ2 P 2 , ( 20
)
where P 2 is the Legendre polynomial of degree two defined by Eq. ( 71). The parameter ᾱ2 is in the range [-1, 2] and is negative (resp. positive) for shells thinner at the equator (resp. poles). The parameterization of the thickness variation is discussed in more detail in Section 4.1.

ACCEPTED MANUSCRIPT

The stress functions for contraction and despinning then satisfy the following duality:

F C w0 = 1 1 -ν -ᾱ2 F D w2 . ( 21
)
The first term in the right-hand side is the contraction solution when the elastic thickness is constant. The duality (21) can be checked by substitution in Eq. ( 14) and using the identity (69). Since spherical harmonics of degree one belong to the kernel of ∆ (see A), the duality remains valid if ᾱ includes such harmonics. If ᾱ contains harmonics up to degree k > 2, a similar relation exists (see Eq. ( 125)) between the solutions for right-hand sides of higher degree ( w = w3 P 3 , ... , w = wk P k ); I use it in F for the computation of the contraction solution at first order. I also give below a relation of fourth degree between the stresses in order to discuss the combination of contraction and despinning.

Eq. ( 21) leads to a duality between the stress resultants for the contraction and despinning solutions:

N C ij w0 = δ ij 1 -ν -ᾱ2 N D ij w2 . ( 22
)
The stress resultants for a contracting shell are thus the same as those for a despinning shell with w2 = -ᾱ2 w0 , plus the stress resultants for a contracting shell of constant thickness. If the contracting shell ( w0 < 0) is thinner at the equator (ᾱ 2 < 0), the equivalent 'despinning' shell is actually spinning faster since w2 < 0. This duality is clearly visible by comparing Figs. 4a and 4b. A similar duality holds for the stresses, but the first term in the right-hand side now depends on latitude:

σC ij w0 = ᾱ δ ij 1 -ν -ᾱ2 σD ij w2 . ( 23
)
Finally the duality relation for the strains reads

C ij w0 = ᾱ δ ij -ᾱ2 D ij w2 . ( 24
)
These dualities are extremely useful. Their most obvious use is to generate both contraction and despinning solutions by computing only one of them, but they also have other advantages. First, Eq. ( 23) shows that the ordering in size of the components of the stress is either the same (if ᾱ2 > 0) or reversed (if ᾱ2 < 0) when going from the contraction to the despinning solution. Thus the azimuthal stress is more compressive than the meridional stress for a despinning shell whether the shell is thinner at the equator or at the poles. While the isotropic term ᾱδ ij /(1ν) in Eq. ( 23) does not affect the ordering of the stress components, its spatial dependence influences the positions of the maxima which must thus be computed from the numerical solution.

Second, the dualities serve to compute a first approximation of the contraction solution on a shell of variable thickness, knowing the despinning solution on a shell of constant thickness. The trick comes from the fact that the term linear in ᾱ2 in the contraction solution is related by duality to the term independent of ᾱ2 in the despinning solution. The contraction solution is computed in such a way in F (see Eqs. ( 121)-( 124)). In dimensional notation, the contraction stresses at first order in ᾱ2 read

σ C θθ ∼ = σ 0 1 + 3 ᾱ2 ν + (2 + ν) cos 2θ 2(5 + ν) , ( 25
)
σ C ϕϕ ∼ = σ 0 1 + 3 ᾱ2 1 + (1 + 2ν) cos 2θ 2(5 + ν) , ( 26
)
where σ 0 is the stress for a radial contraction of δR (δR > 0) when the thickness is constant:

σ 0 = - E 1 -ν δR R . ( 27
)
Fig. 13a shows that the first order approximation is good if ᾱ2 = -2/5 and rather bad if ᾱ2 is close to its lower bound -1. In F, the method is generalized to an ᾱ having an arbitrary Legendre expansion (see Eqs. ( 128)-( 129)). Though the basic characteristics of stress and strain for a despinning shell could be predicted by the symmetry arguments of Section 3.1, duality cannot serve for the computation of their numerical values. Their formulas at first order in ᾱ2 are given in J.

A third interesting consequence of the dualities concerns the analysis of combined contraction and despinning. If the amounts of contraction and despinning are related by w2 = ᾱ2 w0 , (28) the total stress computed with the duality ( 23) is isotropic (though not homogeneous):

σC ij + σD ij = ᾱ w0 δ ij 1 -ν . ( 29
)
Consider a contracting ( w0 < 0) and despinning ( w2 > 0) body with a lithosphere thinner at the equator (ᾱ 2 < 0). If there is less (resp. more) despinning than the threshold (28), the meridional stress is more (resp. less) compressive than the azimuthal stress. This fact is useful when determining the orientation of faults for contraction plus despinning. The threshold is reached for a contraction δR

(δR > 0) of δR = h T 2 ᾱ2 (Ω 2 f -Ω 2 i )R 2 3g , ( 30
)
where Ω i (resp. Ω f ) is the initial (resp. final) angular rate, g is the surface gravity and h T 2 is the degree-two displacement Love number. The above equation results from the substitution of Eqs. ( 89) and (92) into Eq. ( 28).

The threshold (28) becomes dependent on latitude when the inverse thickness expansion includes higher harmonics. If ᾱ = 1+ ᾱ2 P 2 + ᾱ4 P 4 , the generalization of the stress duality reads

σC ij w0 + ᾱ2 σD ij w2 + ᾱ4 σE ij w4 = ᾱ δ ij 1 -ν , ( 31
)
where the index E denotes the solution of the membrane equation ( 14) for w = w4 P 4 , which is given at zeroth order in (ᾱ 2 ,ᾱ 4 ) by Eqs. ( 94)-( 95). Since there is actually no deformation of degree four, the third term in the left-hand side of Eq. ( 31) should not be interpreted as a physical stress but rather as a deviation from the duality (23). At the threshold (28), the total stress σT ij = σC ij + σD ij is thus not isotropic anymore:

σT θθ -σT ϕϕ ∼ = ᾱ4 w0 19 + ν (12P 4 -10P 2 -2) . ( 32
)
If ᾱ4 w0 > 0, the meridional stress is more (resp. less) compressive than the azimuthal stress for high (resp. low) latitudes, with the turning point being given by a latitude of arcsin(1/ √ 7) ∼ = 22.2 • (this value is only valid if ᾱ2 , ᾱ4 1).

If ᾱ4 w0 < 0, the high and low latitudes zones are exchanged. In Section 4, I propose a parameterization of the inverse thickness which has an unlimited expansion in Legendre polynomials (see Eq. ( 36)). Fig. 9 then shows the latitudinal dependence of the threshold between north-south and east-west thrust faults for different types of thickness variations. The case described by Eq. ( 32) corresponds to a value of the parameter k close to zero.

Solving the membrane equation

Symmetry and duality arguments have served us well for the determination of the main characteristics of the stress and strain curves, as well as for a first numerical approximation of the contraction solution. It remains however necessary to solve the full membrane equation in order to compute both the contraction and despinning solutions at an arbitrary degree of precision. I will present two methods which are valid for an arbitrary deformation of the surface.

The first one is numerical and has the advantage of minimizing the amount of programming. The second one is semi-analytical, in the sense that the solution can be expressed as a perturbation expansion in Legendre polynomials about the solution for constant shell thickness. It has the advantage of producing explicit formulas for the stresses in which the influence of parameters describing thickness variations clearly appears. Before tackling the methods of resolution, I briefly discuss the parameterization of the thickness variation.

Thickness variation

Since my aim is to predict tectonic patterns without assuming a specific lithospheric structure, I will work with a simple parameterization of the variation of the shell thickness. I choose to parameterize ᾱ instead of the thickness h because it is more convenient for the semi-analytical method. A few assumptions constrain the possible form of the thickness variation. Variations caused by solar insolation and tidal heating are of very large wavelength, so that only slow-varying functions of the colatitude are permissible. Axial symmetry means that the thickness only depends on the colatitude θ, but it also imposes that the derivative of the thickness vanishes at the poles. If there is equatorial symmetry, the derivative of the thickness vanishes at the equator (this restriction can be lifted without changing the methods of resolution). A simple parameterization of the thickness compatible with these constraints is given by Eq. ( 20):

ᾱ(θ) = 1 + ᾱ2 P 2 (cos θ) , ( 33
)
where the dependence on the colatitude θ has been made explicit. I define the equator-to-pole thickness ratio r by

r = h E h P , ( 34
)
where h E is the equatorial thickness and h P the polar thickness. The coefficient ᾱ2 is then given by

ᾱ2 = 2 r -1 r + 2 . ( 35
)
The thickness is positive everywhere if -1 ≤ ᾱ2 ≤ 2. Negative and positive values of ᾱ2 describe elastic shells that are thinner at the equator and at the poles, respectively. Fig. 2 shows various profiles of ᾱ and h/h P when the lithosphere is thinner at the equator: the values r = (1, 1/2, 1/4, 1/10) correspond to ᾱ2 = (0, -2/5, -2/3, -6/7).

[Figure 2 about here.]

Although the parameterization ( 33) is sufficient for the determination of the most important features of tectonic patterns, it suffers from a weakness revealed by Fig. 2b. When the equator-to-pole thickness ratio decreases, the zone of equatorial thinning expands toward the pole. In other words, the parameter ᾱ2 not only affects the equator-to-pole thickness ratio, but it also modifies the relative extension of the thin and thick zones. Therefore, the parameterization (33) is not appropriate for the analysis of the influence of the thin zone size. I thus define a new parameterization in which the interval [0, π/2] is non-linearly stretched by a function ψ k . As a result, the extension of the thin zone can be modulated without affecting the equator-to-pole thickness ratio:

α(θ) = 1 + α2 P 2 (cos ψ k (θ)) , ( 36
)
where -1 ≤ α2 ≤ 2. The 'tilde' notation indicates that the normalization differs from the one used for ᾱ because the Legendre coefficient of degree zero differs from one. The term P 2 (cos ψ k (θ)) has generally an unlimited expansion in Legendre polynomials, including a term of degree zero. Nonetheless the normalization of α does not matter for the computation of stress and strain; it only affects the stress resultants.

The function

ψ k (θ) is defined on the interval [0, π/2] by ψ k (θ) = ⎧ ⎨ ⎩ π 2 sinh(2kθ/π) sinh(k) if k > 0 , π 2 tanh(2kθ/π) tanh(k) if k < 0 . (37)
It can extended to the interval [0, π] by mirror symmetry about the equatorial plane. The limit k → 0 yields the parameterization (33):

lim k→0 ψ k (θ) = θ . ( 38
)
Fig. 3 shows that positive values of k correspond to a reduction in the size of the thin zone. The equator-to-pole thickness ratio is left unchanged since it is controlled by α2 . Conversely, negative values of k correspond to an extension of the thin zone.

[Figure 3 about here.]

Numerical method

The software Mathematica [START_REF] Research | Mathematica, Version[END_REF] has a powerful command called 'NDSolve' for the numerical solution of differential equations. When there is axial symmetry, the second membrane equation ( 14) is an ordinary differential equation of order four. I found it convenient to solve it as a system of two differential equations of order two:

∆ F = -q , ( 39
) -∆ (ᾱ q) -(1 + ν) A ᾱ ; F = ∆ w . ( 40
)
where the nondimensional load q is defined by

q = R α 0 q . ( 41
)
The various terms in the above equations can be expanded with the formulas of A:

∆ F = ∆ F + 2 F , (42) ∆ (ᾱ q) = (∆ᾱ) q + 2 ∂ ᾱ ∂θ ∂ q ∂θ + ᾱ (∆q) + 2ᾱ q , ( 43
)
A ᾱ ; F = ∂ 2 ᾱ ∂θ 2 + ᾱ cot θ ∂ F ∂θ + F + cot θ ∂ ᾱ ∂θ + ᾱ ∂ 2 F ∂θ 2 + F , ( 44
) with ∆f = ∂ 2 f ∂θ 2 + cot θ ∂f ∂θ . ( 45
)
In principle the system (39)-( 40) should be solved on the segment [0, π] but the points 0 and π must be excluded from the interval because of the apparent polar singularities in spherical coordinates. I thus restrict the interval to [a, b] with a = and b = π -( is a small positive number, for example 0.001).

Four boundary conditions are required. Smooth solutions with axial symmetry satisfy

∂ F ∂θ | θ=a = ∂ q ∂θ | θ=a = 0 , ( 46
) ∂ F ∂θ | θ=b = ∂ q ∂θ | θ=b = 0 . ( 47
)
Before applying these equations as boundary conditions, I should make a caveat.

The stress function has the particularity that it can be redefined at will because of the degree-one freedom mentioned in A:

F → F + β cos θ , ( 48
)
where β is an arbitrary real number. Since cos θ satisfies Eqs. ( 46)-( 47) in the limit → 0, another condition must be specified in order to fully determine F , for example:

F | θ=a = 0 . ( 49
)
This last condition being arbitrary, one should keep in mind that the same physical problem can be solved with various boundary conditions. The resulting stress functions only differ by their degree-one content, i.e. by cos θ multiplied by some constant. The corresponding stresses are of course identical.

Suitable boundary conditions for the interval [a, b] consist of Eq. ( 49), plus three conditions chosen among Eqs. ( 46)-( 47). I found it best to specify three conditions on F and one condition on q. In most cases NDSolve also works well with only Eqs. ( 46)-( 47) as boundary conditions, meaning that its algorithm chooses one possible solution for F among an infinity (note that the freedom in the definition of F only appears in the limit → 0). Though it is not necessary, one can get rid of the cos θ component of F by orthogonalization:

F → F - 3 2 π 0 F cos θ sin θ dθ . ( 50
)
If there is equatorial symmetry, it is sufficient to consider the interval [a, b] with a = and b = π/2. In that case, Eqs. ( 46)-(47) form a suitable set of boundary conditions since cos θ violates Eq. ( 47) when b = π/2. The cos θ component of F then vanishes.

Semi-analytical method

If the elastic thickness is constant, the membrane equation ( 14) is diagonal in the basis of spherical harmonics. The contraction solution is then of harmonic degree zero whereas the despinning solution is of degree two (see C). This straightforward method cannot be used when the elastic thickness is spatially variable because of the coupling of the different harmonic degrees and orders.

In this section I show that the system of coupled differential equations remains manageable under the assumption of axial symmetry. Since the imposed deformation and the elastic thickness do not depend on the longitude, only zonal spherical harmonics contribute.

The problem thus consists in solving the membrane equation ( 14 146)-( 147)). If the expansions are truncated at degree n, the membrane equation takes a matrix form with each row corresponding to a harmonic degree on which the membrane equation is projected, except for degree one:

M (n) F = ∆ w , ( 51
)
with the membrane matrix M (n) being defined by

M (n) = C (n) -(1 + ν) A (n) . ( 52
)
C (n) and A (n) are matrices approximating at degree n the operators C and A, respectively. F and w are vectors containing the coefficients F and w , respectively. ∆ is the diagonal matrix with elements δ = -( + 1) + 2.

The absence of a row for degree one is due to the fact that spherical surface harmonics of degree one do not belong to the images of the operators C and A (see H). Correspondingly, there is no coefficient of degree one in the vectors F and w since the membrane equation does not constrain them (see A). Given a deformation w, the membrane equation in its matrix form can for example be solved by matrix inversion, yielding the Legendre coefficients of the stress function:

F = M (n) -1 ∆ w , ( 53
)
An example of the membrane matrix and of its solution is given in I for an expansion of ᾱ limited to degree two and a truncation degree equal to 6. If the thickness variation is symmetric about the equatorial plane, only Legendre poly-nomials of even degree will contribute to the solution because the deformation is also symmetric (degree zero or two).

Since ᾱ = 1+ ᾱ2 P 2 is finite even if h E or h P vanishes (though not both), the nondimensional membrane equation can in principle be solved with a vanishing thickness at the equator or at the pole. Problems of divergence however appear in these extreme cases for the nondimensional stress function and the stress resultants, though the stresses and the strains are well behaved. The membrane matrix M (n) is invertible in the physical range of ν ∈ [0, 1/2] and for ᾱ2 ∈ [-1, 2],

as can be seen from the examination of its eigenvalues.

If ᾱ is limited to degree two, the approximate membrane equation ( 51) can also be solved by expanding F in ᾱ2 and solving order by order:

F = ∞ p=0 (ᾱ 2 ) p F(p) . ( 54
)
It is convenient to split the membrane matrix so that the dependence in ᾱ2

becomes explicit:

M (n) = M (n) 0 + ᾱ2 M (n) 2 . ( 55
)
The matrix M

(n) 0

is diagonal so that its inversion is straightforward. The solution is then generated order by order in ᾱ2 by a recurrence relation:

F(p+1) = -ᾱ2 M (n) 0 -1 M (n) 2 F(p) , ( 56
)
which is initiated with F(0) , the solution of M

(n) 0 F = ∆ w. This method can be generalized to an ᾱ of degree higher than two. Since the largest eigenvalues of the iteration matrix appearing in Eq. ( 56) tend to 1 (from below) as the truncation degree increases, the series (54) may diverge for near extremal values of ᾱ2 . If |ᾱ 2 | < 1, the series converges but the convergence can be very slow if ᾱ2 is close to -1. Even if it is not the best numerical method, the ᾱi -expansion provides an explicit solution which is helpful to understand properties such as the rule governing the decrease of the Legendre coefficients of the solution (see I) or the pseudo-nodes in the stress and strain components (see F and J).

ACCEPTED MANUSCRIPT

Tectonic patterns

Stress and faulting

Tectonic patterns can be predicted from the analysis of stresses. According to Anderson's theory of faulting, the faulting style depends on how the vertical (or radial) stress compares with the horizontal (or tangential) stresses [START_REF] Melosh | Global tectonics of a despun planet[END_REF]Jaeger et al., 2007;[START_REF] Schultz | Fault populations[END_REF]. Thrust faults, strike-slip faults and normal faults occur if the radial stress is respectively least compressive, intermediate or most compressive among the principal stresses. Thrust faults and normal faults strike in the direction of the intermediate principal stress, whereas strike-slip faults strike in a direction at about 30 • from the direction of the most compressive stress. The radial stress vanishes since faults occur at the surface (this assumption is criticized by [START_REF] Golombek | Fault type predictions from stress distributions on planetary surfaces -Importance of fault initiation depth[END_REF]). Anderson's theory presupposes that all principal stresses are compressive, at least when including the lithostatic pressure that must be added to the stresses computed from thin shell theory. I will briefly mention below the possible occurrence of near-surface tensile failure when one (or more) principal stress is extensional. This point is discussed in more depth by [START_REF] Melosh | Global tectonics of a despun planet[END_REF] and [START_REF] Schultz | Observations, models, and mechanisms of failure of surface rocks surrounding planetary surface loads[END_REF].

Once the contraction solution of the membrane equation (2) has been found with the methods of Section 4, it is possible to compute stress resultants, stresses and strains with the formulas of Section 2. The nondimensional stress is sufficient for the determination of tectonic patterns since the physical magnitude of the stress is only useful when comparing the predicted stresses with a failure criterion. Nonetheless figures represent dimensional quantities for values of parameters typical of terrestrial planets: E = 100 GPa, ν = 0.25, average inverse thickness equal to 1/(100 km), contraction of 0.1% and flattening reduction of 0.15%. I initially give examples with ᾱ limited to harmonic degree two, i.e. ᾱ is parameterized by Eq. ( 33). Most basic characteristics of the tectonic patterns are indeed independent of the precise latitudinal variation of the elastic thickness as long as it smoothly varies from the pole to the equator. Figs. 4 and5 show the stress resultants and the stresses for the contraction and despinning of a lithosphere thinner at the equator. The duality between contraction and despinning solutions appears clearly on Fig. 4. The stress curves for varying ᾱ2 seem to cross at a fixed point on Fig. 5 but this is only true at first order in ᾱ2 (see F and J). The case of polar thinning is illustrated by Fig. 1, which shows the strains for expansion and despinning (stress curves are very similar).

Besides I will resort to the parameterization (36) in order to demonstrate two effects due to the variation in size of the thin zone: it modifies the position of the extremum of the meridional stress and it affects the orientation of thrust faults when contraction is combined with despinning. Fig. 6 shows the stresses for contraction and despinning when the thin zone size is extended or reduced.

[Figure 4

Contraction or expansion

For the contraction of a shell thinner at the equator, the stresses have the following characteristics:

• stresses are everywhere compressive.

• the meridional stress is more compressive than the azimuthal stress (Fig. 5a); This situation favors the development of thrust faults striking east-west.

• the stress is most compressive at the equator (Fig. 5a). This situation favors the development of thrust faults at the equator. Note that the stress resultant is most compressive at the poles (Fig. 4a).

The contraction fault pattern predicted with Anderson's theory is thus significantly modified by a variable elastic thickness: instead of thrust faults with random strikes, the predicted pattern consists of thrust faults striking east-west and preferably formed at the equator. The third characteristic is not strictly true for all thickness profiles. The curve k = -2 in Fig. 6a shows that two things happen if the thin zone has a large extension. First stress curves become nearly flat in that zone, where the stresses approach the value predicted for a lithosphere of constant thickness. Second the latitude of the most compressive meridional stress moves away from the equator toward the latitude where the thick zone begins. This peak in compressive stress is however not wellpronounced and for all practical purposes the most compressive stress can be said to occur in the whole thin zone. If the thin zone is very localized, stress curves become flat in the thick zone and approach the value predicted for a lithosphere of constant thickness (this is not yet apparent on Fig. 6a because k should be much larger than 2).

If the shell is expanding (see Fig. 1a), stresses change sign, so that Anderson's theory predicts normal faults striking east-west and preferably formed at the equator. Near the surface tensile failure is also possible because lithostatic pressure might be too small to render all stresses compressive. It would then lead to the appearance of joints striking east-west, because rocks are much weaker in extension than in compression.

If the contracting shell is thinner at the poles, the azimuthal stress is more compressive than the meridional stress and they are both most compressive at the poles, so that the tectonic pattern consists of thrust faults striking northsouth and preferably forming near the poles. If an expanding shell is thinner at the poles (see Fig. 1a), the azimuthal stress is more extensional than the meridional stress and they are both most extensional at the poles, so that the tectonic pattern consists of normal faults striking north-south and preferably forming near the poles.

Despinning

For the despinning of a shell thinner at the equator, the stresses have the following characteristics:

• the meridional stress is everywhere extensional, whereas the azimuthal stress is extensional near the poles and compressive near the equator. This situation favors the development of normal faults near the poles and strikeslip faults near the equator. The limit between the faulting provinces is given by the latitudes for which the azimuthal stress vanishes. These boundaries are at about ±48 • latitudes (independent of ν) when the elastic thickness is constant, and move toward the poles if the elastic thickness decreases at the equator.

• the meridional stress is more extensional -or less compressive -than the azimuthal stress. This situation either favors normal faults striking eastwest, or strike-slip faults striking at about 60 • from the north (depending on whether the tangential stresses have the same sign or not).

• the meridional stress is most extensional at the poles if the elastic thickness is constant. The maximum moves toward the equator if the elastic thickness variation is large enough (if ᾱ is limited to degree two, the maximum moves away from the pole if ᾱ2 < -0.44). If the equator-to-pole thickness ratio is small (threshold ≈ 1/4) and the thin zone is not too large (threshold k ≈ 1 if equator-to-pole thickness ratio=1/10), the most extensional meridional stress can be at the equator but it is a rather extreme situation (see Fig. 6b).

The despinning fault pattern predicted with Anderson's theory is thus not substantially altered by the variation of the elastic thickness. There is an equatorial province of strike-slip faults striking at about 60 • from the north, plus two polar provinces of normal faults striking east-west. The boundaries between the faulting provinces move toward the poles as the elastic thickness becomes thinner at the equator. As in the case of expansion, tensile failure may occur near the surface, leading to the production of east-west joints. In the rather extreme case where the meridional stress is at the equator, these joints would first form at the equator.

If the shell is spinning up, stresses are the same except for a sign change, so that the tectonic pattern consists of an equatorial province of strike-slip faults striking at about 30 • from the north and polar provinces of thrust faults striking east-west.

If the shell is thinner at the poles (see Fig. 1b), the tectonic pattern is similar to the one for a constant elastic thickness, except that the boundaries between the tectonic provinces shift toward the equator as the shell becomes thinner at the poles.

Contraction plus despinning

Let us now consider tectonics due to simultaneous contraction and despinning. I initially assume that the lithosphere is thinner at the equator and I use the parameterization (33) to start with. Results for the case of polar thinning are given at the end. Fig. 7 (inspired by Fig. 5 of [START_REF] Melosh | Global tectonics of a despun planet[END_REF]) represents the latitudes of the boundaries between tectonic provinces when a varying amount of contraction (or expansion) is added to despinning (Poisson's ratio is chosen to be 0.25 but only weakly affects the position of the boundaries). The proportion between contraction and despinning is parameterized by the contraction/despinning ratio χ:

χ = - w0 w2 , ( 57
)
which is zero if there is only despinning, positive if there is additional contraction, negative if there is additional expansion. Fig. 8 illustrates some possible tectonic patterns. The predicted tectonic pattern has the following features: preferred latitude for the boundary between these provinces and then read on the x-axis the required contraction/despinning ratio χ. Using the same symbols as in Eq. ( 30), I can relate the contraction to the rotation variation:

•
δR = -χ h T 2 (Ω 2 f -Ω 2 i )R 2 3g . ( 58
)
The above equation results from the substitution of Eqs. ( 89) and (92) into

Eq. ( 57). The important conclusion is that a particular combination of con- [Figure 7 about here.]

[Figure 8 about here.]

[Figure 9 about here.]

Applications

Iapetus' ridge

Cassini images revealed in December 2004 an extraordinary feature on Saturn's satellite Iapetus, a high ridge running along the equator spanning more than half of its circumference [START_REF] Porco | Cassini Imaging Science: Initial Results on Phoebe and Iapetus[END_REF]. High-resolution imaging is not available on the whole surface, but it seems that impact basins are present where the ridge has not been detected [START_REF] Denk | Iapetus Imaging During the Targeted Flyby of the Cassini Spacecraft[END_REF]. It is thus reasonable to assume that the ridge originally spanned the whole circumference. In wellimaged areas, it has an average width of 60 -70 km and heights up to 18 km [START_REF] Giese | The topography of Iapetus' leading side[END_REF][START_REF] Denk | Iapetus Imaging During the Targeted Flyby of the Cassini Spacecraft[END_REF].

Several theories have been proposed to explain the origin of the ridge, including deposition of ring remnants, convection, despinning and compaction.

The morphologic analysis of the ridge by [START_REF] Giese | The topography of Iapetus' leading side[END_REF] excluded an exogenous origin due to the accretion of a now-vanished ring [START_REF] Ip | On a ring origin of the equatorial ridge of Iapetus[END_REF]. Convection [START_REF] Czechowski | The Iapetus's ridge: Possible explanations of its origin[END_REF] how this surface reduction could be localized along the equator.

The tectonic patterns studied in this paper provide new possibilities. In particular, the contraction of a lithosphere thinner at the equator leads to a tectonic pattern of compressional faults striking east-west and preferably formed at the equator. The two assumptions underlying this mechanism are not arbitrary.

First, the lithosphere must be thinner at the equator because of the latitudinal variation in solar insolation, as shown in Figs. Though it is not known whether faulting occurred at the ridge, other tectonic processes such as buckling or folding will result in the same orientation and location since the meridional stress is more compressive than the azimuthal stress and the maximum compression is at the equator.

How does despinning affect the contraction hypothesis? Iapetus is at present locked in a 1:1 resonance with Saturn: its rotation period and orbital period are both equal to 79 days. This synchronicity is generic of all large satellites which initially rotated much faster. Moreover Iapetus has the peculiarity of being very flattened (equatorial and polar radii differ by 35 km), though it should be nearly spherical if it were hydrostatic. Its present shape is well fitted by a homogeneous hydrostatic ellipsoid with a rotation period of about 16 hours [START_REF] Castillo-Rogez | Iapetus' geophysics: Rotation rate, shape, and equatorial ridge[END_REF]. The most obvious explanation is that Iapetus had an initial rotation period of 16 hours or less and that its shape froze while despinning. The lithosphere must have been several hundred kilometer thick when the rotation period reached 16 hours, otherwise the body would not have retained its flattened shape. In the evolution model of [START_REF] Castillo-Rogez | Iapetus' geophysics: Rotation rate, shape, and equatorial ridge[END_REF], despinning is a rapid process occurring 800 Myr after formation, at which time the lithosphere has a thickness of about 240 km. In this scenario, contraction happens long before despinning, when the lithosphere is thin enough to be easily deformed. In such a case the east-west faulting pattern due to contraction is not directly affected by despinning.

If contraction occurred before despinning, a question to be addressed is the absence of a superimposed tectonic pattern due to despinning. This problem is not unique to Iapetus. Though all large satellites underwent despinning, none exhibits unambiguous evidence of it in its global tectonic grid. As for Iapetus, one way out is to postulate an initial period not much shorter than 16 hours. If the lithosphere is thick enough, the flattening change is small and the associated stresses not large enough to cause faulting.

How thick should the lithosphere be in order to resist despinning stresses?

The simplest failure criterion is that of Coulomb (Jaeger et al., 2007), which consists of a linear relation between the extreme principal stresses. For thrust and normal faults, this criterion can be rewritten as a constraint between the maximum differential stress σ f = |σ 1σ 3 | (σ 1 and σ 3 are the most and least compressive stresses, respectively) and the vertical stress σ z = ρgz:

σ f = a + b σ z , ( 59
)
where a is the yield strength of ice at zero pressure and σ z is positive. If σ z ≤ 12 MPa, the constants (a, b) for ice at a temperature of about 100 K are equal to (3.4 MPa, 1.86) and (1.2 MPa, 0.65) for compression and extension, respectively [START_REF] Beeman | Friction of ice[END_REF]. Failure starts either at the poles or at the equator since the differential stress is maximum at these locations. At the poles, normal faults start to form when the meridional stress exceeds the yield strength of ice for extension. At the equator, thrust faults form instead of strikeslip faults when the shell is thick (see D): faulting begins when the azimuthal stress exceeds (in absolute value) the yield strength of ice for compression. As discussed in Section 2, a thick lithosphere can usually be approximated by a shell of constant thickness so that Eqs. ( 108 7 and8, the combination of contraction and despinning leads to three possible faulting patterns at the equator, given in order of increasing contraction: strike-slip faults, north-south thrust faults, or east-west thrust faults. The amount of contraction required to yield the third outcome is given by Eq. ( 30), at least if the lithospheric thickness varies as in Eq. ( 33). As a starting point, I set the initial rotation period to 16 hours, the equator-to-pole thickness ratio to 1/2 (i.e. ᾱ2 = -2/5) and the Love number h T 2 to 5/2 (limit of vanishing lithospheric thickness). The contraction threshold between northsouth and east-west thrust faults is then δR ≈ 60 km which is large, though of the same order of magnitude as the one due to porosity collapse in the model of [START_REF] Castillo-Rogez | Iapetus' geophysics: Rotation rate, shape, and equatorial ridge[END_REF]. Of course, the amount of contraction due to porosity depends on the initial porosity profile which is unknown. Other values of the parameters strongly affect the result. On the one hand, the threshold increases with the square of the initial rotation rate. On the other hand the threshold decreases as the lithospheric thickness increases, because the flattening variation is then smaller. The two-layer incompressible model discussed in D yields δR ≈ 13 km and δR ≈ 0.9 km if the lithospheric thickness is equal to 20 km and 240 km, respectively. Furthermore the threshold not only decreases with the equator-to-pole thickness ratio (see Fig. 7b), but is also affected by the precise form of the thickness variation: Fig. 9b shows that a reduction of the thin zone can decrease the value of the threshold at the equator by a factor of two. Therefore, simultaneous contraction and despinning may lead, or not, to east-west thrust faults. For a given amount of contraction, the outcome sensitively depends on the initial rotation rate and on the lithospheric thickness (average thickness, equator-to-pole thickness ratio, thin zone size).

Besides the contraction hypothesis, an episode of expansion in Iapetus' evolution due to differentiation [START_REF] Squyres | The tectonics of icy satellites[END_REF] would lead to equatorial normal faults striking east-west and forming first at the equator. Iapetus' ridge does not look like normal faulting, but the same extensional stress field is compatible with the intrusive dike origin for the ridge proposed by [START_REF] Melosh | An Intrusive Dike Origin for Iapetus' Enigmatic Ridge?[END_REF]. In that scenario, the simultaneous occurrence of despinning is generally not a problem because it is only required that the meridional stress be most extensional at the equator. This is true for moderate amounts of expansion, the value of the threshold depending on the thickness profile. The expansion threshold is for example 7.5 times smaller in magnitude than the contraction threshold discussed above, assuming an equator-to-pole thickness ratio of 1/2 (it coincides in this particular case with the expansion threshold for normal faults at the equator). Note that despinning can also be favorable to the formation of an east-west dike at the equator, but the conditions are rather stringent: the equator-to-pole thickness ratio must be smaller than 1/4 and the thin zone not too large (see Fig. 6b).

Mercury's scarps

Lobate scarps on Mercury are linear or arcuate escarpments varying in length from about 20 km to 500 km and in height from a few hundred meters to about 3 km [START_REF] Strom | Tectonism and volcanism on Mercury[END_REF]Watters et al., 2009). They are the most common tectonic features on Mercury and are thought to be the result of thrust faulting.

The predominant north-south orientation of lobate scarps was first attributed to illumination bias but [START_REF] Watters | Thrust faults and the global contraction of Mercury[END_REF] showed that this preferred orientation is real. The combination of planetary contraction and despinning is a possible explanation for this pattern [START_REF] Melosh | Mercurian global tectonics -A consequence of tidal despinning[END_REF][START_REF] Pechmann | Global fracture patterns of a despun planet -Application to Mercury[END_REF][START_REF] Dombard | Despinning plus global contraction and the orientation of lobate scarps on Mercury: Predictions for MESSENGER[END_REF]. Another interesting observation is that lobate scarps in the southern polar region are as often east-west than north-south (see Fig. 2.23 in Watters and Nimmo (2009)). East-west lobate scarps are much rarer in the northern polar region, but the total length of scarps of any orientation is also much smaller than in the south. Let us suppose that lobate scarps predominantly strike north-south from 50 • S to 50 • N latitude and east-west in the polar regions, as done by [START_REF] King | Pattern of lobate scarps on Mercury's surface reproduced by a model of mantle convection[END_REF]. Fig. 9a shows that simultaneous contraction and despinning may lead to this pattern if the thin zone is rather extended. As mentioned in the introduction, the lithosphere could be thinner by a factor of two at the equator in comparison with the poles [START_REF] Mckinnon | Application of ring tectonic theory to Mercury and other solar system bodies[END_REF][START_REF] Melosh | The tectonics of Mercury[END_REF]. In that case, the required contraction/despinning ratio should be close to two, assuming that the shape parameter k is between -1 and -2. Substituting χ = 2 in Eq. ( 58), I obtain δR ≈ 20 km if h T 2 = 5/2 and the initial rotation period is 20 hours (g = 3.7 m/s 2 , R = 2440 km). The contraction threshold is somewhat reduced by more realistic values of h T 2 . The thickness of the lithosphere at the time of the formation of lobate scarps is difficult to pinpoint: arguments have been made for values of 25 -30 km and 100 -200 km (see a summary in [START_REF] Watters | The Tectonics of Mercury[END_REF]). In any case, Fig. 10a shows that the secular Love number h T 2 is not very sensitive to the lithospheric thickness. The contraction threshold is only reduced to 15 km and 12 km if the lithospheric thickness is equal to 30 km and 100 km, respectively (see D for the value of h T 2 ). The radius contraction at the origin of lobate scarps likely being in the range 1-2 km [START_REF] Strom | Tectonism and volcanism on Mercury[END_REF]Watters et al., 2009), it is far too small to make possible the existence of east-west thrust faults.

One kilometer of contraction, combined with despinning, is also too small to generate north-south thrust faults all over the planet, though two kilometers are enough if the yield strength of rock is neglected (use Eq. ( 58) with χ determined from Fig. 7a). Yet another problem is that lobate scarps are more abundant in the south than near the equator, contrary to the prediction of contraction plus despinning. [START_REF] Dombard | Despinning plus global contraction and the orientation of lobate scarps on Mercury: Predictions for MESSENGER[END_REF] thus suggested that a larger amount of contraction occurred during despinning in the early history of Mercury (their contraction threshold of 3-5.5 km is however significantly reduced if the factor h T 2 is taken into account). The Late Heavy Bombardment then erased the faults and new scarps were produced along the old fault lines by a later global contraction event of smaller magnitude, when the planet had despun. Such a scenario widens the choice of the relative weights of despinning and contraction.

On the one hand despinning needs not be finished when the early tectonic grid forms. On the other hand early contraction can be rather large in thermal evolution models: 2.5 -7.5 km in [START_REF] Dombard | Despinning plus global contraction and the orientation of lobate scarps on Mercury: Predictions for MESSENGER[END_REF], with an upper bound of 17 km if the core completely solidifies [START_REF] Solomon | Some aspects of core formation in Mercury[END_REF]) (measurements of Mercury's librations however indicate that the core is at least partially molten, see [START_REF] Margot | Large Longitude Libration of Mercury Reveals a Molten Core[END_REF]). In conclusion, it is possible to imagine a sequence of events in which various events of contraction and despinning lead to thrust faults that strike north-south near the equator and east-west near the poles.

However the observation that the total scarp length increases from north to south remains to be explained. East-west thrust faults could instead be the result of the reactivation of early normal faults due to despinning [START_REF] Watters | The Tectonics of Mercury[END_REF]. The various orientations of lobate scarps could also be due to a combination of despinning, contraction and true polar wander [START_REF] Matsuyama | Gravity and tectonic patterns of Mercury: Effect of tidal deformation, spin-orbit resonance, nonzero eccentricity, despinning, and reorientation[END_REF].

Conclusions

The main result of the paper can be very simply stated: the contraction or expansion of a planet with lithospheric thinning at the equator results in tectonic features striking east-west, preferably formed in the equatorial region. If the lithosphere is thinner at the poles, contraction or expansion generates tectonic features striking north-south. In general, despinning alone cannot produce eastwest features near the equator. One exception could be the appearance of eastwest joints or dikes at the equator if the thickness variation is large and strongly peaked at the equator.

A spectacular and rather unique illustration of an east-west structure is the equatorial ridge on Iapetus, which can be explained by the theory developed in this paper. The mountainous structure of the ridge is suggestive of a compressional event but an extensional process is not excluded. Since despinning inhibits the formation of east-west compressional tectonics, the ridge probably formed long before Iapetus had despun, but this is not required if the contraction is large enough.

Contraction and despinning stresses can be computed either numerically or with a semi-analytical method, which yields a perturbation expansion in the parameters describing the thickness variation. First order formulas for stresses are sufficient if the thickness varies by a factor of two. In the simplest case (inverse thickness of degree two), contraction stresses are well approximated by the elementary formulas ( 25)-( 26), while Eqs. ( 128)-( 129) are valid at first order for an arbitrary thickness profile.

The combination of radius change and despinning generates new tectonic patterns when the lithospheric thickness is variable. On the one hand, contraction plus despinning makes possible the coexistence of tectonic provinces of thrust faults differing in orientation if the lithosphere is thinner at the equator.

Lobate scarps on Mercury roughly follow such a pattern but the small amount of contraction associated with the scarps makes it necessary to resort to rather complicated models involving fault reactivation. On the other hand, expansion plus despinning may lead to the coexistence of tectonic provinces of normal faults differing in orientation if the lithosphere is thinner at the poles. In both cases the magnitude of the contraction/despinning ratio must be large otherwise the predicted patterns are similar to those valid for a lithosphere of constant thickness.

Global contraction or expansion occur in many models of interior evolution and are often assumed to be the cause of planetary tectonics. Since lithospheric thinning due to the latitudinal variation in solar insolation must be a generic phenomenon, it is surprising that an east-west orientation of tectonic features at the planetary scale is so rare in the solar system. One reason could be that the thinning is too small to have a significant effect. Also, tidal heating might counteract the effect of the variation in solar insolation. Another explanation is the simultaneous occurrence of despinning. Finally, it is possible that global contraction or expansion events generally happen very early in the history of the planet so that either faults do not appear because the lithosphere is not yet formed, or faults do form but the associated tectonic pattern is subsequently erased. Further inferences about the origin of tectonic features depend on more complete tectonic mapping, which is underway for Mercury with Messenger data and for Saturn's icy satellites with Cassini data. efited from exchanges with Surdas Mohit, Francis Nimmo and James Roberts. I also thank Attilio Rivoldini and Antony Trinh for discussions on Love numbers.

A. Differential operators on the sphere

The operators O i are linear differential operators of the second degree on the sphere:

O 1 = ∂ 2 ∂θ 2 + 1 , ( 60
)
O 2 = csc 2 θ ∂ 2 ∂ϕ 2 + cot θ ∂ ∂θ + 1 , ( 61
)
O 3 = cscθ ∂ 2 ∂θ∂ϕ -cot θ ∂ ∂ϕ , ( 62
)
where θ is the colatitude and ϕ is the longitude. They give zero when applied on spherical surface harmonics of degree one.

The scalar differential operator ∆ is defined by [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]:

∆ = O 1 + O 2 (63) = ∆ + 2, ( 64
)
where ∆ is the spherical Laplacian (called surface Laplacian in [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]).

Spherical surface harmonics of degree are eigenfunctions of ∆ with eigenvalues (Blakely, 1995, pp. 121-122)

δ = -( + 1) + 2 . ( 65
)
The scalar differential operators C and A appearing in the membrane equation (2) are defined by

C(a ; b) = ∆ (a ∆ b) , (66) A(a ; b) = (O 1 a)(O 2 b) + (O 2 a)(O 1 b) -2 (O 3 a)(O 3 b) . ( 67
C (1 ; b) -(1 + ν) A (1 ; b) = (∆ -1 -ν) ∆ b . (68) C (a ; 1) -(1 + ν) A (a ; 1) = (1 -ν) ∆ a . ( 69
)

B. Operators O i on Legendre polynomials

Legendre polynomials P (x) are defined in many books [START_REF] Whittaker | A Course of Modern Analysis[END_REF][START_REF] Blakely | Potential Theory in Gravity and Magnetic Applications[END_REF]. I only give those that appear in the text, in the form P (cos θ):

P 0 = 1, ( 70
)
P 2 = 1 4 (3 cos 2θ + 1) , ( 71
)
P 4 = 1 64 (35 cos 4θ + 20 cos 2θ + 9) , ( 72
)
P 6 = 1 512 (231 cos 6θ + 126 cos 4θ + 105 cos 2θ + 50) . ( 73
)
The action of the operators O i defined by Eqs. ( 60)-( 62) on Legendre polynomials P (cos θ) is

O 1 P = 1 -2 P + Σ , ( 74
)
O 2 P = (1 -) P -Σ , ( 75
)
O 3 P = 0, (76)
where

Σ = p r=1 (2 -4 r + 1) P -2r , = (2 -3) P -2 + (2 -7)P -4 + ... ( 77
)
with p = /2 for even and p = ( -1)/2 for odd . The action of O 2 is obtained from an expression for the first derivative of a Legendre polynomial found in Whittaker and Watson (1935, p. 330). The action of O 1 is then obtained from Eqs. ( 63) and ( 65).

The values of O 1 P and O 2 P are equal at the poles:

O i P | θ=0 = 1 2 δ (i = 1, 2) , ( 78
)
O i P | θ=π = (-1) 2 δ (i = 1, 2) , ( 79
)
where δ is defined by Eq. ( 65).

Since the coefficients multiplying the Legendre polynomials in Eqs. ( 74)-( 75) are independent of θ, the operators O i can be represented as matrices acting on vectors of Legendre coefficients. The action of O 1 on a vector of Legendre coefficients of degrees 0, 2, 4 and 6 is represented by the matrix

O (6) 1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 1 1 1 0 -3 5 5 0 0 -15 9 0 0 0 -35 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 80
)
The action of O 2 on a vector of Legendre coefficients of degrees 0, 2, 4 and 6 is represented by the matrix

O (6) 2 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 -1 -1 -1 0 -1 -5 -5 0 0 -3 -9 0 0 0 -5 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 81
)

C. Constant thickness: thin shell

When the shell thickness h is constant, the membrane equation ( 2) becomes

(∆ -1 -ν) ∆ F = Eh R ∆ w , ( 82
)
where I used the identity (68). Eq. ( 82) is diagonalized by expanding F and w in spherical surface harmonics whose coefficients ( = 1) satisfy (see Eq. ( 65))

F m = Eh R 1 -( + 1) + 1 -ν w m , (83)
The index m refers to the order which is equal to zero under the assumption of axial symmetry. The component of degree one remains indeterminate. The stress function generates the stresses according to Eqs. ( 4)-( 5). The action of O i on Legendre polynomials of degree 0, 2, 4, 6 is given by Eqs. ( 80)-( 81).

For a contracting shell, w = w 0 P 0 and the stress function is constant (P 0 = 1):

F C = Eh 1 -ν w 0 R P 0 . ( 84
)
The non-zero stresses are

σ C θθ = σ C ϕϕ = E 1 -ν w 0 R . ( 85
)
For a despinning shell, w = w 2 P 2 and the stress function reads

F D = - Eh 5 + ν w 2 R P 2 . ( 86
)
The non-zero stresses are

σ D θθ = - E 5 + ν w 2 R (-P 2 -1) , ( 87
)
σ D ϕϕ = - E 5 + ν w 2 R (-3P 2 + 1) . ( 88
)
The coefficient w 2 can be related to the flattening variation δf (δf < 0 for despinning) by

w 2 = - 2 3 R δf . ( 89
)
After substitution of Eqs. ( 71) and ( 89), despinning stresses [START_REF] Vening-Meinesz | Shear patterns of the Earth's crust[END_REF][START_REF] Melosh | Global tectonics of a despun planet[END_REF] read

σ D θθ = - E 5 + ν δf 6 (3 cos 2θ + 5) , ( 90
) σ D ϕϕ = - E 5 + ν δf 6 (9 cos 2θ -1) . ( 91
)
The flattening is related to the angular rate Ω by

f = h T 2 Ω 2 R 2g , ( 92
)
where g is the surface gravity and h T 2 is the degree-two tidal Love number for radial displacement [START_REF] Matsuyama | Tectonic patterns on reoriented and despun planetary bodies[END_REF].

For a deformation of degree four, w = w 4 P 4 and the stress function reads

F E = - Eh 19 + ν w 4 R P 4 . ( 93
)
The non-zero stresses are

σ E θθ = - E 19 + ν w 4 R (-3P 4 -5P 2 -1) , ( 94
)
σ E ϕϕ = - E 19 + ν w 4 R (-15P 4 + 5P 2 + 1) . ( 95
)

D. Constant thickness: thick shell

In a tidal or centrifugal potential, the radial deformation of the surface of a spherically stratified planet is parameterized by the Love number h T 2 , as in Eq. ( 92). Other Love numbers noted l T 2 and k T 2 parameterize the tangential deformation of the surface and the gravity perturbation at the surface induced by a tidal or centrifugal potential [START_REF] Lambeck | The earth's variable rotation: Geophysical causes and consequences[END_REF]. Love numbers are complicated functions of the internal planetary structure. A simple formula exists for an incompressible body of uniform density and shear modulus (see below), but this model is not appropriate for long-term processes such as despinning. Stresses indeed relax by viscous flow or creep everywhere within the planet except in the lithosphere, which is by definition rigid on geologic time scales. Internal layers in which stresses have relaxed are modeled as fluid, inertial forces are neglected (static limit) and the resulting Love numbers are called secular [START_REF] Matsuyama | Tectonic patterns on reoriented and despun planetary bodies[END_REF]. The simplest interior model including a lithosphere is a two-layer spherically stratified planet made of an incompressible fluid core surrounded by an incompressible elastic lithosphere [START_REF] Melosh | Global tectonics of a despun planet[END_REF]. A more realistic density profile is obtained by dividing the planet into several homogeneous layers, each of constant thickness. Tidal Love numbers for such models are readily computed with the propagator matrix technique [START_REF] Sabadini | Global Dynamics of the Earth[END_REF]. The assumption of incompressibility is not realistic but satisfactory for our purpose, since the main uncertainty regarding secular Love numbers lies in the unknown thickness of the lithosphere.

In the limit of vanishing core radius, the two-layer model gives the degreetwo tidal Love numbers for a homogeneous incompressible body (Love, 1944, pp. 257-259):

h T 2 = 5 2 1 + 19 2 µ ρgR -1 , ( 96
) l T 2 = 3 10 h T 2 , ( 97
)
k T 2 = 3 5 h T 2 . ( 98
)
R is the radius, ρ is the density, g is the surface gravity and µ is the shear modulus. Note that [START_REF] Munk | The rotation of the earth; a geophysical discussion[END_REF] give an incorrect formula for l T 2 while Lambeck (1980) gets it right. Love numbers for a hydrostatic homogeneous planet are obtained by setting µ = 0: h T 2 = 5/2 and k T 2 = 3/2. The Love number l T 2 does not exist for a hydrostatic body since displacements in a fluid are indeterminate in the static limit, apart from the radial displacement at the surface.

In the limit of vanishing elastic thickness, the two-layer incompressible model gives the following degree-two tidal Love numbers:

h T 2 = 5 2 , ( 99
)
l T 2 = 3 11 h T 2 , ( 100
) k T 2 = 3 5 h T 2 . ( 101
)
As expected, h T 2 and k T 2 coincide with the values for a homogeneous hydrostatic body. More generally, a lithosphere of vanishing thickness neither affects the shape of the body nor its gravitational field, whatever the internal density distribution. However l T 2 differs from the value for a homogeneous incompressible body since the presence of a membrane modifies tangential displacements at the surface.

I now show that the value of l T 2 given by Eq. ( 100) agrees with the membrane limit of thin shell theory. Toroidal displacement is absent when the shell thickness is constant (axial symmetry is actually sufficient), assuming that there is no toroidal external potential. The tangential displacement v can then be expressed as the surface gradient of a potential S [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]:

v = ∇ S .
(102)

The displacement potential S is related to the stress function by Eq. ( 71) of [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF], which becomes in the membrane limit

∆ S = Rα (1 -ν)∆ F -2w . ( 103
)
The substitution in this equation of the stress function ( 83) yields the degreecomponent of S in terms of w :

S = 1 + ν ( + 1) -1 + ν w . ( 104
) If V 2 = -Ω 2 R 2 /3
is the degree-two component of the external potential at the surface, S 2 and w 2 are related to degree-two tidal Love numbers by

w 2 = h T 2 g V 2 , ( 105
)
S 2 = l T 2 g V 2 .
(106)

In the membrane limit, thin shell theory thus predicts that

l T 2 = 1 + ν 5 + ν h T 2 . ( 107
)
This relation holds for any density distribution within the planet, as long as it is spherically symmetric. If the lithosphere is incompressible (ν = 1/2), the ratio l T 2 /h T 2 is equal to 3/11, that is Eq. ( 100). If the lithosphere is made of ice (ν ≈ 1/3), the same ratio is equal to 1/4 (this result has been noted by [START_REF] Wahr | Modeling stresses on satellites due to nonsynchronous rotation and orbital eccentricity using gravitational potential theory[END_REF] regarding tidal stresses).

The insertion of Eq. ( 89) and Eqs. ( 105)-( 106) into stress-strain relations yields the following expressions for the tangential stresses at the surface of a spherically stratified planet:

σ D θθ = - E 1 -ν 2 2 3 δf (1 + ν) P 2 + l T 2 h T 2 1 -ν -2(2 + ν)P 2 , ( 108
)
σ D ϕϕ = - E 1 -ν 2 2 3 δf (1 + ν) P 2 + l T 2 h T 2 -1 + ν -2(1 + 2ν)P 2 , ( 109
)
where I used the fact that the rr component of the stress vanishes at the surface.

In the membrane limit of thin shell theory, Love numbers satisfy Eq. ( 107) so that Eqs. ( 108)-( 109) become equivalent to the stresses ( 90)-( 91). Moreover the substitution of Eq. ( 97) into Eqs. ( 108)-( 109) yields the same stresses as those found by [START_REF] Melosh | Global tectonics of a despun planet[END_REF] for a homogeneous incompressible body (see his

Eq. ( 20)).

Once Love numbers are known for some interior model, it becomes possible to compute surface stresses and thus predict despinning tectonic patterns using

Anderson's theory of faulting (see Section 5). This is a convenient way to compute thick shell effects on faulting style as was first done by [START_REF] Melosh | Global tectonics of a despun planet[END_REF], since several codes for the computation of Love numbers are now available. The boundaries of tectonic provinces depend only on two numbers: the ratio l T 2 /h T 2 and Poisson's ratio ν. For example, despinning induces an equatorial zone of north-south thrust faults if the meridional stress is compressional at the equator, that is if

l T 2 h T 2 < 1 + ν 6 . ( 110
)
This inequality is satisfied neither by thin shell theory (see Eq. ( 107)) nor by the homogeneous incompressible model (see Eq. ( 97)), but is verified for the incompressible models shown on Fig. 10b if the shell thickness is between 8.4% and 62% of the radius. Though h T 2 and l T 2 are very sensitive to the internal structure (see Fig. 10a), their ratio varies much less. That is why the despinning tectonic pattern is not much affected by the shell thickness.

[Figure 10 about here.] For Iapetus, I adopt an incompressible two-layer model with a radius of 735.6 km and a uniform density of 1083 kg/m 3 in both layers [START_REF] Thomas | Shapes of the saturnian icy satellites and their significance[END_REF]. The shear modulus µ of ice at 257 K is 3.52 GPa so that its value at 80 K (approximately the mean temperature of Iapetus' surface) is about 4.4 GPa [START_REF] Gammon | Elastic constants of artificial and natural ice samples by Brillouin spectroscopy[END_REF]). Poisson's ratio ν for ice is 0.325 (but ν = 0.5 in the model) and Young's modulus E = 2(1 + ν)µ is 11.7 GPa. The model yields h T 2 = 0.53 and h T 2 = 0.038 when the lithospheric thickness is equal to 20 km and 240 km, respectively (for other values, see Fig. 10).

For Mercury, I adopt an incompressible three-layer model with a fluid core, fluid mantle and elastic lithosphere. Parameters are chosen to be in the middle range of admissible values [START_REF] Rivoldini | The interior structure of Mercury and its core sulfur content[END_REF]: core radius of 2000 km, density of mantle and lithosphere equal to 3100 kg/m 3 (thus core density equal to 7325 kg/m 3 ), seismic wave velocities in the lithosphere given by v p = 7900 m/s and v s = 4550 m/s. The corresponding elastic parameters for the lithosphere are: µ = 64 GPa, E = 160 GPa, ν = 0.25 (but ν = 0.5 in the model). The model yields h T 2 = 1.83 and h T 2 = 1.49 when the lithospheric thickness is equal to 30 km and 100 km, respectively (for other values, see Fig. 10).

E. Two properties of the strain

The tangential displacement in the azimuthal direction vanishes if there is axial symmetry, i.e. if the elastic thickness and the radial deformation of the shell only depend on the latitude. In that case, the strain-displacement equations (e.g. [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]

read θθ = 1 R ∂v θ ∂θ + w , ( 111
) ϕϕ = 1 R (cot θ v θ + w) , ( 112
) θϕ = 0, ( 113
)
where w is the radial displacement and v θ the tangential displacement in the meridional direction.

By symmetry, v θ vanishes at the poles. The average of θθ on [0, π] is thus related to the average of w = w/R:

< θθ > = < w > , ( 114
)
where

< f >= 1 π π 0 f dθ . ( 115
)
There is equatorial symmetry if the elastic thickness and the radial deformation are symmetric about the equatorial plane. In that case, v θ vanishes at the equator so that the average can be done on [0, π/2]. For contraction plus despinning, w = w0 + w2 P 2 , so that

< θθ > = w0 + w2 4 . (property 1) (116)
Another consequence of axial symmetry is that the azimuthal strain becomes independent of the elastic thickness at the equator:

ϕϕ | θ= π 2 = w| θ= π 2 (117)
This property also directly results from the variation in length of the equator, δL = 2π w| θ= π 2 , which is uniformly distributed along the equator. For contraction plus despinning,

ϕϕ | θ= π 2 = w0 - w2 2 .
(property 2) (118)

F. Contraction solution at first order

Given the despinning solution on a shell of constant thickness, the dualities of Section 3.2 generate a first approximation of the contraction solution on a shell of variable thickness. If ᾱ depends only on the parameter ᾱ2 as in Eq. ( 20), the contraction and despinning solutions can be expanded in ᾱ2 about the solution for a shell of constant thickness given in C. The dualities ( 21)-( 24) then allow to compute the contraction solution at order ᾱn+1 2 in terms of the despinning solution at order ᾱn 2 . The nondimensional despinning stress function at zeroth order is obtained from Eq. ( 86):

F D w2 = - P 2 5 + ν + O(ᾱ 2 ) . ( 119
)
Nondimensionalization has been carried out with Eqs. ( 9)-( 13). The duality (21) yields the contraction stress function at first order in ᾱ2 :

F C w0 = 1 1 -ν + ᾱ2 P 2 5 + ν + O(ᾱ 2 2 ) . ( 120
)
The duality (23) acting on Eqs. ( 87)-( 88) yields the contraction stresses at first order in ᾱ2 (an alternative is to use Eqs. ( 15)-( 16) on F C ):

σC θθ w0 ∼ = 1 1 -ν + ᾱ2 2(2 + ν)P 2 -1 + ν (1 -ν)(5 + ν) , ( 121
)
σC ϕϕ w0 ∼ = 1 1 -ν + ᾱ2 2(1 + 2ν)P 2 + 1 -ν (1 -ν)(5 + ν) . ( 122
)
Fig. 13a compares Eqs. ( 121)-( 122) with the stresses computed at high order in ᾱ2 , that is with a truncation degree n = 20. The first order approximation is good for small ᾱ2 (ᾱ 2 = -2/5) but is bad for large ᾱ2 (ᾱ 2 = -6/7).

The strains at first order in ᾱ2 can be computed with Eqs. ( 17)-( 19):

C θθ w0 ∼ = 1 + ᾱ2 1 + ν 5 + ν (4P 2 -1) , ( 123
) C ϕϕ w0 ∼ = 1 + ᾱ2 1 + ν 5 + ν (2P 2 + 1) . ( 124
)
With the values P 2 | θ=0 = 1, P 2 | θ=π/2 = -1/2 and < P 2 >= 1/4, one can check that the stresses and strains ( 121)-( 124) satisfy the properties enumerated at the beginning of Section 3.1. Besides the meridional strain as given by Eq. ( 123) is equal to its average value at θ = 45 • (where P 2 = 1/4), so that the meridional strain curves for varying ᾱ2 seem to have a unique crossing point (see Fig. 1a).

Since this property is only valid at first order in ᾱ2 , I call this point a pseudonode. Pseudo-nodes for stresses are obtained by setting to zero the term of order ᾱ2 in Eqs. ( 121)-( 122), yielding respectively θ ∼ = 48.2 • and θ ∼ = 65.9 • (ν = 0.25)

for the meridional and azimuthal stresses (see Fig. 5a).

The same method works for ᾱ not limited to degree two. The generalization of the duality ( 21) to an ᾱ depending on an arbitrary number of parameters ᾱ reads

F C w0 = 1 1 -ν - ≥2 ᾱ F L w , ( 125
)
where the index L denotes the solution of the membrane equation ( 14) for w = w P . At first order in ᾱ , the contraction stress function can thus be computed from the solutions (83) on a shell of constant thickness:

F C w0 ∼ = 1 1 -ν + ≥2 ᾱ P f ν , ( 126
)
where

f ν = ( + 1) -1 + ν . ( 127
)
There is no contribution of ᾱ1 to the contraction stress function ( 126), but this is only valid at first order in ᾱ (terms in ᾱn 1 vanish, but mixed terms like ᾱ1 ᾱ2 do not). At first order in ᾱ , the contraction stresses read

σC θθ w0 ∼ = 1 1 -ν + ≥1 ᾱ O 2 P f ν + P 1 -ν , ( 128
) σC ϕϕ w0 ∼ = 1 1 -ν + ≥1 ᾱ O 1 P f ν + P 1 -ν , ( 129
)
with the action of O i on P defined by Eqs. ( 74)-( 75). The sum now starts at = 1 because of the multiplication by ᾱ when computing the stresses from the stress resultants, but the first term in the parentheses vanishes for = 1. At first order in ᾱ , the contraction strains read

C θθ w0 ∼ = 1 -(1 + ν) ≥1 ᾱ f ν (O 1 -1) P , ( 130
) C ϕϕ w0 ∼ = 1 -(1 + ν) ≥1 ᾱ f ν (O 2 -1) P . ( 131
)
The stresses ( 128)-( 129) and the strains ( 130)-( 131) satisfy the properties enumerated at the beginning of Section 3.1 because of Eqs. ( 78)-( 79) and since

Eqs. ( 60)-( 61) with axial symmetry imply

< (O 1 -1) P > = 0, ( 132
) (O 2 -1) P | θ= π 2 = 0. (133)

G. Operators C and A on spherical surface harmonics

If the shell thickness is constant, the membrane equation reduces to Eq. ( 82)

which depends only on the operator ∆ and can be solved with spherical surface harmonics as in C. If the thickness is variable, the membrane equation is not diagonal in the basis of spherical surface harmonics. This basis remains however useful if the action of the operators C and A on spherical surface harmonics generates a linear combination of the same functions, so that the operators can be represented as square matrices. The action of the operators O i -from which C and A are built -on Legendre polynomials is known (see Eqs. ( 74)-( 76))

and could be used to derive a formula for the action of C and A on Legendre polynomials. I will however obtain a more compact formula for the action of A by rewriting it exclusively in terms of the scalar operator ∆ . The operator A can be written in the equivalent form [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]:

A(a ; b) = (∆a)(∆ b) -(∇ i ∇ j a)(∇ i ∇ j b) + (∆ a) b + a (∆ b) + 2 a b , ( 134
)
where ∇ i is the covariant derivative on the sphere and summation on repeated indices is implicit (indices are raised with the inverse metric g ij = diag(1, sin -2 θ)).

Repeated covariant differentiation yields the following identities:

∆(ab) = (∆a) b + 2 (∇ i a)(∇ i b) + a (∆ b) (135) ∆∆(ab) = (∆∆a) b + 2 (∆ a)(∆ b) + a (∆∆ b) + 4 (∇ i ∇ j a)(∇ i ∇ j b) + 4 (∇ i a)(∇ i b) + 4 (∇ i ∆ a)(∇ i b) + 4 (∇ i a)(∇ i ∆b) . ( 136
)
In the derivation of the last identity, I used the commutation relation for the spherical Laplacian and the covariant derivative on the surface of the sphere,

∆∇ i a -∇ i ∆ a = ∇ i a , ( 137
)
which can be derived from the general commutation relation of covariant derivatives of a vector (see Eq. (G1) in [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]. Other useful identities are obtained by substituting a → ∆ a or b → ∆ b in Eq. ( 135). With the above identities, the operator A can be rewritten as an expression where the only intervening operator is ∆ :

A(a ; b) = 1 4 [ -∆ ∆ (ab) -(∆ ∆ a) b -a (∆ ∆ b) + 2 (∆ a)(∆ b) + 2 ∆ ((∆ a) b + a (∆ b)) -2 ∆ (ab) -2 (∆ a) b -2 a (∆ b) + 8 ab ] . ( 138
)
The action of the operators C and A can be partially evaluated if a and b are spherical surface harmonics of degree m and n (noted a m and b n ):

C(a m ; b n ) = δ n ∆ (a m b n ) , ( 139
)
A(a m ; b n ) = - 1 4 ( ∆ ∆ + κ mn ∆ + λ mn ) (a m b n ) . ( 140
)
where δ n is given by Eq. ( 65) and

κ mn = 2(1 -δ m -δ n ) , ( 141
)
λ mn = (δ m -δ n ) 2 + 2 (δ m + δ n ) -8 . ( 142
)
The action of the operators C and A can only be fully evaluated if the product a m b n is expanded in spherical surface harmonics.

H. Operators C and A on Legendre polynomials

The action of the operators C and A on spherical surface harmonics a m and b n is given by Eqs. ( 139)-( 140) in which products a m b n appear. I am now going to evaluate these products with the assumption that the spherical surface harmonics are zonal or, equivalently, Legendre polynomials.

The product of two Legendre polynomials can be expanded into Legendre polynomials with the Adams-Neumann formula (Whittaker and Watson, 1935, p. 331):

P m (x) P n (x) = min(m,n) j=0 k mnj P m+n-2j (x) , ( 143
)
with the coefficients of the expansion being given by

k mnj = m + n + 1 2 -2j m + n + 1 2 -j K j K m-j K n-j K m+n-j , ( 144
)
where

K m = (2m -1)!! m! = 1.3...(2m -1) m! , ( 145
)
with K 0 = 1. These coefficients k mnj are symmetric in (m, n) and k m00 =

k 0m0 = 1.
The action of the operators C and A on Legendre polynomials of degree m and n is thus (the dependence on x is omitted):

C(P m ; P n ) = min(m,n) j=0 k mnj χ mnj P m+n-2j , ( 146
)
A(P m ; P n ) = min(m,n) j=0 k mnj ψ mnj P m+n-2j , (147)
where

χ mnj = δ n δ m+n-2j , ( 148
)
ψ mnj = - 1 4 δ m+n-2j 2 + κ mn δ m+n-2j + λ mn , ( 149
)
with δ n given by Eq. ( 65) and coefficients (κ mn , λ mn ) given by Eqs. ( 141)-( 142).

If m = 0 (constant elastic thickness), the only relevant coefficients are

χ 0n0 = (δ n ) 2 , ( 150
)
ψ 0n0 = δ n . ( 151
)
The Legendre polynomial of degree one does not belong to the images of the operators C and A: if m + n -2j = 1 (meaning that |m -n| = 1), then δ m+n-2j = 0 and λ mn = 0, so that χ mnj = ψ mnj = 0. The more general proof that spherical surface harmonics of degree one are not included in the images of C and A is given in [START_REF] Beuthe | Thin elastic shells with variable thickness for lithospheric flexure of oneplate planets[END_REF]. Since the coefficients multiplying the Legendre polynomials in Eqs. ( 146)-( 147) are independent of θ, the operators C and A can be represented as matrices acting on vectors of Legendre coefficients.

I. Membrane matrix for n = 6

The membrane matrix formalism described in Section 4 is illustrated here with a simple example: ᾱ is limited to degree two and the truncation degree n is equal to 6. The action of the operator C on the vector F = F0 , F2 , F4 , F6 T can be computed with Eq. ( 146). The matrix C (6) appearing in Eq. ( 52) is given by the sum

C (6) = C (6) 0 + ᾱ2 C (6)
2 , where

C (6) 0 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 4 k 000 0 0 0 0 16k 020 0 0 0 0 324 k 040 0 0 0 0 1600 k 060 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , C (6) 2
= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 -8 k 222 0 0 -8 k 200 16 k 221 72 k 242 0 0 72k 220 324 k 241 720 k 262 0 0 720 k 240 1600 k 261 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
The action of the operator A on the vector F is computed with Eq. ( 147). The matrix A (6) appearing in Eq. ( 52) is given by the sum

A (6) = A (6) 0 + ᾱ2 A (6)
2 , where

A (6) 0 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 2 k 000 0 0 0 0 -4 k 020 0 0 0 0 -18 k 040 0 0 0 0 -40 k 060 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , A (6) 2
= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 -4 k 222 0 0 -4 k 200 20 k 221 6 k 242 0 0 6k 220 90 k 241 24 k 262 0 0 2 4k 240 200 k 261 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .
The four rows are associated with the projection of the membrane equation on harmonic degrees 0, 2, 4, 6. The upper index on the matrices denotes that the expansion is limited to degree 6. The integer coefficients are generated by the functions χ mnj and ψ mnj defined by Eqs. ( 148)-( 149). In particular, the integer coefficients in the matrices C (6) 0

and A (6) 0 are given by Eqs. ( 150)-( 151). The Adams-Neumann coefficients k mnj , given by Eq. ( 144), are in symbolic form so that their origin is apparent: the mn indices indicate that the term comes from the coefficients of ᾱ and F of degree m and n, respectively, whereas m+n-2j is equal to the harmonic degree associated with the row of the matrix. because it is associated with a polynomial of degree 8.

The approximate membrane equation ( 51) now takes the following form:

M (6) F0 , F2 , F4 , F6 T = (2 w0 , -4 w2 , 0, 0) T . ( 155
)
The coefficients F C of the contraction solution up to ᾱ3 2 are generated with the recurrence relation (56), starting with F(0) = ( w0 /(1ν), 0, 0, 0):

F C 0 w0 ∼ = 1 1 -ν + 2 5 1 5 + ν ᾱ2 2 + 4 35 1 + 5ν (5 + ν) 2 ᾱ3 2 , ( 156
) F C 2 w0 ∼ = 1 5 + ν ᾱ2 + 2 7 1 + 5ν (5 + ν) 2 ᾱ2 2 + 12 7 49 -7ν + 17ν 2 + ν 3 (5 + ν) 3 (19 + ν) ᾱ3 2 , ( 157
) F C 4 w0 ∼ = 6 35 ν -11 (5 + ν)(19 + ν) ᾱ2 2 + 36 385 (ν -11)(-21 + 56ν + 5ν 2 ) (5 + ν) 2 (19 + ν) 2 ᾱ3 2 , ( 158
) F C 6 w0 ∼ = 18 385 (ν -11)(ν -29) (5 + ν)(19 + ν)(41 + ν) ᾱ3 2 . ( 159
)
Coefficients F C have a leading term of order ᾱ /2

2 . Thus the truncation at harmonic degree n can be understood as an expansion in ᾱ2 : harmonic coefficients that are neglected are at most of order ᾱ1+n/2 2 : if n = 6, F C with ≥ 8 are at most of order O(ᾱ 4 2 ) and are ignored. However the truncation at degree n is not equivalent to an expansion at order ᾱn/2 2 since the lower the degree of the harmonic coefficient, the better it is known. More precisely, the band-diagonal structure of the membrane matrix has the effect that, for a given truncation degree n, the expansion of a term F C can be done up to order ᾱn-/2+1 2 included:

if n = 6, ( F C 0 , F C 2 , F C 4 , F C 6 ) can be expanded up to order (ᾱ 7 2 , ᾱ6 2 , ᾱ5 2 , ᾱ4 2 ) respec- tively.
The coefficients F D of the despinning solution are related to the coefficients of the contraction solution by Eq. ( 21):

F D 0 w2 = - 1 ᾱ2 F C 0 w0 - 1 1 -ν , ( 160
)
F D i w2 = - 1 ᾱ2 F C i w0 (i = 2, 4, 6) . ( 161
)
Coefficients F D thus have a leading term of order ᾱ /2-1 2

, except F D 0 which begins at order ᾱ2 . Furthermore, for a given truncation degree n, coefficients Since the stresses and the strains include a supplementary factor ᾱ vanishing at the poles when ᾱ2 = -1, they are not affected by the convergence problem.

Poisson's ratio has generally a weak effect on the solutions for contraction and despinning, except for the coefficients = 0 for contraction and to a lesser extent = 2 for despinning. The other coefficients are mainly affected by ν at near extremal values of ᾱ2 .

[Figure 11 about here.]

[Figure 12 about here.]

The computation of the stresses completes the example. The nondimensional stress resultants are computed with Eq. ( 15) and Eqs. ( 74)-( 76):

Nθθ = F0 + F2 O 2 P 2 + F4 O 2 P 4 + F6 O 2 P 6 , ( 163
) Nϕϕ = F0 + F2 O 1 P 2 + F4 O 1 P 4 + F6 O 1 P 6 . ( 164
)
The Legendre coefficients of O i P j are given by the columns of the matrices (80)-(81). Nθϕ is zero. The nondimensional stress is related to Nij by Eq. ( 16):

σij = (1 + ᾱ2 P 2 ) Nij . ( 165
)
The product of Legendre polynomials can be computed with the Adams-Neumann formula (143) but I will not give the explicit formulas. The stresses and the strains for the contraction solution at first order in ᾱ2 have already been given in F whereas those for the despinning solution are computed in J.

J. Despinning solution at first order

If ᾱ depends only on the parameter ᾱ2 , the despinning stress function at first order in ᾱ2 can be obtained from the contraction solution ( 156)-( 158) and the dualities ( 160)-( 161):

F D w2 ∼ = - P 2 5 + ν - ᾱ2 5 + ν 2 5 + 2 7 1 + 5ν 5 + ν P 2 + 6 35 ν -11 19 + ν P 4 .
The despinning stresses are computed at first order in ᾱ2 from F D with Eqs. ( 15)-( 16):

σD θθ w2 ∼ = P 2 + 1 5 + ν + ᾱ2 s θθ (ν, P 2 ) (5 + ν) 2 (19 + ν) . ( 166
) σD ϕϕ w2 ∼ = 3P 2 -1 5 + ν + ᾱ2 s ϕϕ (ν, P 2 ) (5 + ν) 2 (19 + ν) . ( 167
)
where the auxiliary functions defined by

s θθ (ν, x) = -31 + 18ν + ν 2 + 69 + 48ν + 3ν 2 x + 40 + 18ν + 2ν 2 x 2 , s ϕϕ (ν, x) = 21 -30ν -3ν 2 + 47 + 72ν + ν 2 x + 10 + 42ν + 8ν 2 x 2 ,
have the property

s θθ (ν, 1) = s ϕϕ (ν, 1) , ( 168
)
which ensures the equality of the stresses at the poles. In Eqs. ( 166)-( 167), P 4 has been expressed in terms of P 2 :

P 4 = 35 18 (P 2 ) 2 - 2 7 P 2 - 1 5 . ( 169
)
The strains can be computed at first order in ᾱ2 with Eqs. ( 17)-( 19):

D θθ w2 ∼ = (1 -3ν)P 2 + 1 + ν 5 + ν + ᾱ2 (1 -ν 2 ) e θθ (ν, P 2 ) (5 + ν) 2 (19 + ν) . ( 170
) D ϕϕ w2 ∼ = (3 -ν)P 2 -1 -ν 5 + ν + ᾱ2 (1 -ν 2 ) e ϕϕ (ν, P 2 ) (5 + ν) 2 (19 + ν) . ( 171
)
where the auxiliary functions defined by

e θθ (ν, x) = -31 -3ν + (69 + ν) x + (40 + 8ν) x 2 , ( 172
) e ϕϕ (ν, x) = (2x + 1) (21 + ν + (5 + ν) x) , (173)
have the properties e θθ (ν, 1) = e ϕϕ (ν, 1) , (174)

< e θθ (ν, P 2 ) > = 0, (175)

e ϕϕ (ν, -1/2) = 0 , ( 176
)
which ensure that the properties enumerated at the beginning of Section 3.1 are satisfied (note that < P 2 >= 1/4 and < (P 2 ) 2 >= 11/32).

At zeroth order in ᾱ2 , the meridional strain is equal to its average value at θ = 45 • (where P 2 = 1/4), but this is not true at first order in ᾱ2 since e θθ (ν, 1/4) = 0. If ν = 0.25, the pseudo-nodes defined in Section F occur at θ ∼ = 40.3 • for the meridional strain (see Fig. 1b), θ = 43.8 • for the meridional stress and θ = 64.4 • for the azimuthal stress (see Fig. 5b).

Fig. 13b compares Eqs. ( 166)-( 167) with the stresses computed at high order in ᾱ2 , that is with truncation degree n = 20. The first order approximation is good for small ᾱ2 (ᾱ 2 = -2/5) but is rather bad for large ᾱ2 (ᾱ 2 = -6/7).

[Figure 13 about here.]

Faulting style and orientation for simultaneous despinning and contraction (or expansion) when the lithosphere is thinner at the equator. The thickness variation is parameterized by the equator-to-pole thickness ratio r as on Elastic parameters are E = 100 GPa and ν = 0.25. The thickness variation is parameterized by the equator-to-pole thickness ratio r as on Fig. 2; the average inverse thickness is equal to 1/(100 km). Truncation degree n is equal to 20. Tensile stress is positive. The duality ( 22) is apparent when comparing contraction and despinning curves. : Faulting style and orientation for simultaneous despinning and contraction (or expansion) when the lithosphere is thinner at the equator. The thickness variation is parameterized by the equator-to-pole thickness ratio r as on Fig. 2. The x-axis is the contraction/despinning ratio χ defined by Eq. ( 57). Poisson's ratio is equal to 0.25. (a) Boundaries of the areas with a given fault style. For given ᾱ2 and χ, the intersections of the curves with a vertical line of abscissa χ give the latitudes of the boundaries between tectonic provinces. (b) Contraction threshold beyond which the orientation of thrust faults changes from north-south to east-west. If ᾱ is limited to degree two as is the case here, the threshold is independent of latitude and given by χ = -1/ ᾱ2 . The threshold for r = 1 (ᾱ 2 = 0) is at infinity because the orientation of thrust faults never switches from north-south to east-west when the elastic thickness is constant. The boundaries of the areas with a given fault style are shown as in (a). 2 /h T 2 . The model is described in D. For Mercury, the lithospheric thickness is limited by the size of the core. On an incompressible planet, despinning leads to thrust faults in the equatorial region when l T 2 /h T 2 < 1/4 (see Eq. ( 110) with ν = 0.5); this threshold is indicated by the horizontal line in graph (b).

  ) with expansions in zonal spherical harmonics, i.e. Legendre polynomials: ᾱ = ᾱ P , F = F P and w = w P . The action of the operators C and A on Legendre polynomials produces a finite sum of Legendre polynomials which is computed in G and H (final results are embodied in Eqs. (

  about here.] [Figure 5 about here.] [Figure 6 about here.]

  If there is no contraction or expansion, strike-slip faults are predicted near the equator and east-west normal faults near the poles. The boundaries of tectonic provinces move by a few degrees with respect to the case of constant elastic thickness.• If the despinning planet also contracts, the strike-slip fault province extends toward the poles. Beyond a first threshold of contraction, an area of north-south thrust faults appears near the equator; this area gets larger if contraction increases, whereas the strike-slip fault province is split in two smaller parts which are displaced toward the poles. Beyond a second threshold of contraction, the provinces of normal and strike-slip faults vanish. Thrust faults are then predicted over the whole planet. The two thresholds have a moderate dependence on the equator-to-pole thickness ratio and a weak dependence on Poisson's ratio.• When contraction exceeds the first threshold, thrust faults striking northsouth start to form. Since thrust faults for pure contraction strike eastwest, the orientation of the faults must change for a large enough contraction. If ᾱ is limited to degree two, Eqs. (28)-(29) show that the orientation of thrust faults switches everywhere from north-south to east-west if w2 = ᾱ2 w0 , that is if χ = -1/ ᾱ2 .• If the despinning planet also expands, the strike-slip faulting province becomes smaller and vanishes for a large enough expansion. East-west normal faults are then predicted all over the planet. If the thickness variation is small, normal faults preferably form in polar areas. Beyond some expansion threshold, normal faults preferably form in the equatorial region, except if the thickness is constant.If ᾱ is not limited to degree two, the orientation of thrust faults may depend on the latitude. With the assumptions ᾱ2 , ᾱ4 1 and other ᾱi = 0 (i > 4), Eq. (32) gives the sign of the stress difference at the threshold w2 = ᾱ2 w0 . If ᾱ4 < 0, the change from north-south to east-west thrust faults first occurs in the polar region. As contraction increases, the frontier between provinces of thrust faults with different orientations shifts from the pole to the equator. Conversely, the change from north-south to east-west thrust faults first occurs in the equatorial region if ᾱ4 > 0 and the frontier moves from the equator to the pole as contraction increases. These two cases are illustrated by Figs. 9a,b where the parameterization (36) is used: negative values of ᾱ4 correspond to negative values of k (enlargement of thin zone), whereas positive values of ᾱ4 correspond to positive values of k (reduction of thin zone). For a given k, one can select a

  traction and despinning leads to tectonic provinces of thrust faults differing in orientation.Finally I consider simultaneous contraction and despinning on a lithosphere thinner at the poles. In contrast with the case of lithospheric thinning at the equator, the tectonic pattern when contraction (resp. expansion) is dominant consists of north-south thrust faults (resp. normal faults) preferably formed near the poles. Hence thrust faults do not change in orientation as the amount of contraction increases but the area where they preferably form is displaced from the equator to the poles. The change in orientation occurs in the expansion regime: normal faults switch from striking east-west to north-south when w2 = ᾱ2 w0 (see Eqs. (28)-(29) with ᾱ2 > 0). Near the threshold, provinces of normal faults with different orientations (north-south or east-west) can coexist.

  13 and 14 of[START_REF] Howett | Thermal inertia and bolometric Bond albedo values for Mimas, Enceladus, Tethys, Dione, Rhea and Iapetus as derived from Cassini/CIRS measurements[END_REF].Second, a collapse of porosity during the first 15 Myr of Iapetus occurs in the interior model of[START_REF] Castillo-Rogez | Iapetus' geophysics: Rotation rate, shape, and equatorial ridge[END_REF], leading to a global contraction of about 40 km during which the average lithospheric thickness is less than 20 km.

  ) Since C(a ; b), A(a ; b) and ∆ b give zero if b is a spherical surface harmonic of degree one, the membrane equation (2) does not constrain the degree-one components of the stress function F and of the transverse displacement w (translation invariance).C and A are linear in a and in b. If a is constant A(a ; b) = a∆ b whereas A(a ; b) = b∆ a if b is constant, so that the following identities hold:

F

  Figs. 11 and 12 show the effect of the truncation degree n on the coefficients of the contraction solution. The effects on the despinning solution are similar since they are related by Eqs. (160)-(161). Since the stresses are the quantities of physical interest, the figures do not show the Legendre coefficients of the stress function F but rather those of the sum of the nondimensional stress resultants: N ≡ Nθθ + Nϕϕ = δ F . (162) Truncation effects can be large close to ᾱ2 = 2, but I only consider negative values of ᾱ2 relevant to equatorial thinning. Fig. 11 shows the effect of the truncation degree and of the expansion in ᾱ2 on the first four Legendre coefficients of the solution. The comparison of truncation degrees n = 6 and n = 20 indicates that that the truncation degree n = 6 already yields a precise solution except when ᾱ2 is close to -1. The solution obtained by direct matrix inversion is preferable to the solution obtained by expansion in ᾱ2 . Fig. 12 shows the decrease in the values of the Legendre coefficients of the contraction solution up

Fig. 2 .Figure 2 :Figure 4 :

 224 Fig.2. The x-axis is the contraction/despinning ratio χ defined by Eq. (57).Poisson's ratio is equal to 0.25. (a) Boundaries of the areas with a given fault style. For given ᾱ2 and χ, the intersections of the curves with a vertical line of abscissa χ give the latitudes of the boundaries between tectonic provinces. (b) Contraction threshold beyond which the orientation of thrust faults changes from north-south to east-west. If ᾱ is limited to degree two as is the case here, the threshold is independent of latitude and given by χ = -1/ ᾱ2 . The threshold for r = 1 (ᾱ 2 = 0) is at infinity because the orientation of thrust faults never switches from north-south to east-west when the elastic thickness is constant. The boundaries of the areas with a given fault style are shown as in (a). . . . . . . . . . . . . . . . . . . . . . . 79 8 Tectonic patterns predicted by Anderson's theory for simultaneous despinning and contraction/expansion when the lithosphere is thinner at the equator. This is the pictorial interpretation of Fig. 7. The y-axis is the latitude. The third pattern from left is for despinning only; expansion and contraction are respectively added to the left and to the right of it. Normal faults, strike-slip faults and thrust faults are respectively represented with dashed, dotted and continuous lines. Faulting preferably occurs where the lines are the thickest. The figure does not show all possible transition states. The second pattern from the left is only realized if the equator-to-pole thickness ratio is close to one. The transition pattern EW-NS is shown for the case of an enlarged thin zone (k < 0). . . . . . . . . . . . . . . . . . . . . . . 80 9 Prediction of thrust fault orientation for simultaneous despinning and contraction when the extension of the thin zone varies as on Fig. 3 (r = 1/2; supplementary curves are drawn for k = ±1). These figures generalize Fig. 7b: the vertical line in both graphs corresponds to the vertical line r = 1/2 in Fig. 7b. Poisson's ratio is equal to 0.25. (a) Thin zone enlargement: latitude of the transition between the polar province of east-west faults and the equatorial province of north-south faults. (b) Thin zone reduction: latitude of the transition between the polar province of north-south faults and the equatorial province of east-west faults. . . . . . . . . . . . . . . . . . . 81 10 Secular degree-two tidal Love numbers for Mercury and Iapetus as a function of the lithospheric thickness:(a) h T 2 , (b) l T 2 /h T 2 .The model is described in D. For Mercury, the lithospheric thickness is limited by the size of the core. On an incompressible planet, despinning leads to thrust faults in the equatorial region when l T 2 /h T 2 < 1/4 (see Eq. (110) with ν = 0.5); this threshold is

Figure 5 :Figure 6 :

 56 Figure 5: Stress as a function of colatitude when the lithosphere is thinner at the equator: (a) contraction, (b) despinning. Continuous curves refer to σ θθ while dashed curves refer to σϕϕ. Other details as in Fig. 4.

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure8: Tectonic patterns predicted by Anderson's theory for simultaneous despinning and contraction/expansion when the lithosphere is thinner at the equator. This is the pictorial interpretation of Fig.7. The y-axis is the latitude. The third pattern from left is for despinning only; expansion and contraction are respectively added to the left and to the right of it. Normal faults, strike-slip faults and thrust faults are respectively represented with dashed, dotted and continuous lines. Faulting preferably occurs where the lines are the thickest. The figure does not show all possible transition states. The second pattern from the left is only realized if the equator-to-pole thickness ratio is close to one. The transition pattern EW-NS is shown for the case of an enlarged thin zone (k < 0).

Figure 12 :Figure 13 :

 1213 Figure12: Decrease with harmonic degree of the Legendre coefficients of the contraction solution ( w0 = -1). The coefficients N (see Eq. (162)) are plotted as a function of the harmonic degree (up to = 20) for different values of ᾱ2 . The symbols overlap when = 0. Truncation degree n is equal to 50.

  Apartfrom the trivial values k m00 = k 0m0 = 1, the values relevant to the example above are associated with the products P 2 P 2 , P 2 P 4 and P 2 P 6 :

	(k 220 , k 221 , k 222 ) =	18 35	,	2 7	,	1 5	,	(152)
	(k 240 , k 241 , k 242 ) =	5 11	,	20 77	,	2 7	,	(153)
	(k 260 , k 261 , k 262 ) =	28 65	,	14 55	,	45 143	.	(154)
	The coefficient k 260 does not appear in C	(6) 2	(6) and A 2

< 1/4 (see Eq. (110) with ν = 0.5); this threshold is indicated by the horizontal line in graph (b). . . . . . . . . . . . . . . . 82

Acknowledgments

This work was financially supported by the Belgian PRODEX program managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office. I thank Isamu Matsuyama for suggesting the inclusion of Love numbers and Jay Melosh for discussions on the numerical solution. I ben-

List of Figures 1

Strain as predicted from symmetry arguments: (a) expansion, (b) despinning. Continuous curves refer to θθ while dashed curves refer to ϕϕ. Three cases are shown: constant thickness, equatorial thinning, polar thinning. The unlabeled middle curves in the despinning case are for constant elastic thickness. The numerical values used to draw the curves are: w0 = 0.001 for expansion, w2 = 0.001 for despinning, r = 1/2 for equatorial thinning, r = 2 for polar thinning (r is defined by Eq. ( 34)), Poisson's ratio equal to 0.25.

Extensional strain is positive. . . . . . . . . . . . . . . . . . . . . . . 73 2

Equatorial thinning of the lithosphere as parameterized by Eqs. ( 33)-( 35):

(a) profiles of ᾱ (proportional to the inverse thickness), (b) thickness profiles normalized by the polar thickness. The four profiles correspond to equatorto-pole thickness ratios r equal to (1, 1/2, 1/4, 1/10). . . . . . . . . . . . 74 3

Equatorial thinning of the lithosphere as parameterized by Eqs. ( 36)-( 38):

(a) profiles of α (proportional to the inverse thickness), (b) thickness profiles normalized by the polar thickness. The equator-to-pole thickness ratio r is equal to 1/2: the curves k = 0 are the same as the curves r = 1/2 in Fig. 2a,b Continuous curves refer to θθ while dashed curves refer to ϕϕ. Three cases are shown: constant thickness, equatorial thinning, polar thinning. The unlabeled middle curves in the despinning case are for constant elastic thickness. The numerical values used to draw the curves are: w0 = 0.001 for expansion, w2 = 0.001 for despinning, r = 1/2 for equatorial thinning, r = 2 for polar thinning (r is defined by Eq. ( 34)), Poisson's ratio equal to 0.25. Extensional strain is positive. 162)) are plotted as a function of ᾱ2 for = 0, 2, 4, 6.

The contraction is w0 = -1 and Poisson's ratio is equal to 0.25. The curve for = 0 has been shifted, its real y-intercept being 2/(ν -1) = -8/3 instead of -0.8. Three different approximations to the solution are plotted: n=6 with expansion limited to ᾱ3 2 (dotted curve), n=6 with no expansion in ᾱ2 (dashed curve), n=20 with no expansion in ᾱ2 (continuous curve).