Neil M Ribe
email: ribe@fast.u-psud.fr

Ralf Hielscher

Olivier Castelnau

An analytical finite-strain parametrization for texture evolution in deforming olivine polycrystals

Keywords: Seismic anisotropy, Dynamics of lithosphere and mantle, Mechanics, theory, and modelling, Plasticity, diffusion, and creep

I N T RO D U C T I O N

The principal minerals composing the Earth's uppermost mantle -olivine and orthopyroxene -have strongly anisotropic elastic properties. If the crystals of these minerals contained in a given volume are orientated randomly, the elastic properties of the aggregate will be isotropic. If however the crystals have a non-random orientation distribution (crystallographic preferred orientation or CPO), the speed of an elastic wave traversing the volume will depend on the direction of propagation. The observable consequences of this seismic anisotropy include birefringence (splitting) of shear waves, azimuthal variation of surface wave speeds and differences in propagation speed between surface waves with different polarizations, that is, Rayleigh and Love waves [START_REF] Montagner | Seismic anisotropy and global geodynamics[END_REF][START_REF] Long | Mantle dynamics and seismic anisotropy[END_REF].

The principal cause of CPO is the progressive deformation experienced by mantle rocks as they participate in the global convective circulation. Under appropriate conditions of stress, temperature and grain size, olivine and orthopyroxene crystals deform via the mechanism of dislocation creep, whereby internal dislocations move through the crystal along a limited number of crystallographic planes and directions (slip systems). Deformation of this type constrains the crystallographic axes to rotate relative to a fixed external reference frame, much as a row of tilted books on a shelf rotates when one pushes down on it. Because crystals with different orientations rotate at different rates, the overall distribution of orientations evolves with time in a way that reflects the geometry of the slip systems, the internal stresses that build up within the aggregate during plastic strain and the character of the imposed deformation.

Because CPO and seismic anisotropy are so directly linked to progressive deformation, observations of seismic anisotropy have the potential to constrain the pattern of convective flow in the mantle. Realizing this potential, however, requires a reliable quantitative model for how CPO depends on the history of the deformation experienced by the material. The central question is thus one of polycrystal mechanics: given a polycrystalline aggregate subject to an imposed macroscopic stress or strain rate, how do the individual crystals that compose it deform and rotate? Three general classes of models have been proposed to answer this question.

The first class comprises so-called full-field models. In these, the polycrystal is treated explicitly as a spatially extended body composed of contiguous grains, and the stress and strain within it are field variables that vary continuously as a function of position. Solutions are obtained using either the finite-element method (e.g. [START_REF] Sarma | Effects of interactions among crystals on the inhomogeneous deformation of polycrystals[END_REF] or a more computationally efficient fast Fourier transform (e.g. [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Lebensohn | N-site modeling of a 3D viscoplastic polycrystal using fast Fourier transform[END_REF]. Full-field models allow the stress and strain to vary both among and within individual grains in a physically realistic way. However, their great computational expense makes them too slow (by many orders of magnitude) for routine use in convection calculations.

This disadvantage is overcome to some extent by mean-field homogenization models, in which the detailed spatial distribution of the grains is ignored and the aggregate is treated as a finite number of coexisting mechanical phases having different orientations and (possibly) material properties. Accordingly, compatibility of stress and strain is not enforced between spatially contiguous grains (as it is in full-field models), but rather between each phase and a homogeneous effective medium defined by the average of all the other phases. While homogenization models sacrifice some of the physical realism of full-field models, they are much more computationally efficient. The best-known variant of the general homogenization scheme is the viscoplastic self-consistent (VPSC) model [START_REF] Molinari | A self-consistent approach of the large deformation polycrystal viscoplasticity[END_REF][START_REF] Lebensohn | A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: Application to zirconium alloys[END_REF], which has been widely used to model the evolution of CPO in deforming olivine-dominated aggregates (e.g. [START_REF] Wenk | Modelling plastic deformation of peridotite with the self-consistent theory[END_REF][START_REF] Tommasi | Forward modeling of the development of seismic anisotropy in the upper mantle[END_REF][START_REF] Li | Subduction-induced mantle flow, finite strain, and seismic anisotropy: Numerical modeling[END_REF]. The current state of the art in homogenization models is represented by the second-order self-consistent (SOSC) model developed by Ponte Castañeda (2002) and Liu & Ponte Castañeda (2004), which is discussed further below. Detailed comparisons of full-field and homogenization models have been performed by Lebensohn et al. (2004aLebensohn et al. ( ,b, 2011) ) and (for olivine) by [START_REF] Castelnau | Micromechanical modeling of the viscoplastic behavior of olivine[END_REF].

A final degree of both physical simplicity and computational efficiency is reached in models of the kinematic class, which are based on a quantitative description of the deformation-induced rate of crystallographic rotation as a function of grain orientation. Examples of such models include [START_REF] Ribe | A theory for plastic deformation and textural evolution of olivine polycrystals[END_REF] and the D-Rex model of [START_REF] Kaminski | A kinematic model for recrystallization and texture development in olivine polycrystals[END_REF] and [START_REF] Kaminski | D-Rex, a program for calculation of seismic anisotropy in the convective upper mantle[END_REF], which has been widely used to predict CPO and seismic anisotropy in upper mantle flow models (e.g. [START_REF] Lassak | Seismic characterization of mantle flow in subduction systems: Can we resolve a hydrated mantle wedge? Earth planet[END_REF][START_REF] Conder | Rapid mantle flow beneath the Tonga volcanic arc[END_REF][START_REF] Becker | Radial seismic anisotropy as a constraint for upper mantle rheology[END_REF][START_REF] Faccenda | Seismic anisotropy around subduction zones: Insights from three-dimensional modeling of upper mantle deformation and SKS splitting calculations[END_REF]. An improved kinematic model that matches the predictions of the SOSC model nearly exactly is the ANPAR model of [START_REF] Goulding | Analytical parametrization of self-consistent polycrystal mechanics: Fast calculation of upper mantle anisotropy[END_REF]. Kinematic models are computationally ≈10 4 times faster than homogenization models. However, because the kinematic approach does not account explicitly for stress compatibility among grains, it cannot be used to predict rheological properties of a deforming aggregate.

Despite the important differences between homogenization and kinematic models, they share a common feature: a discrete representation of the aggregate as a finite number N of individual mechanical phases or grains. Typically, large values of N (∼10 3 -10 4 ) are required for statistically acceptable predictions. Moreover, because the orientation of each grain in the aggregate is described by three Eulerian angles, the total number of variables required to describe the texture is 3N. This makes calculations of evolving CPO very cumbersome and time-consuming, even if a relatively fast kinematic model is used; and the problem becomes especially severe if the calculation involves tracking a large number of material particles in a time-dependent flow field. More fundamentally, using 10 3 -10 4 parameters to describe an anisotropic aggregate is superfluous to the point of absurdity given the limited power of seismological data to constrain mantle anisotropy in situ. Shear-wave splitting measurements, for example, can under most conditions constrain only two depth-averaged parameters, the fast propagation direction and the percentage of anisotropy [START_REF] Silver | Shear wave splitting and subcontinental mantle deformation[END_REF].

Motivated by the shortcomings of the traditional methods, we propose here a new approach to calculating texture evolution that combines the physical rigor of homogenization methods with a computational efficiency exceeding that of kinematic models by a factor of ∼10 3 -10 4 . The essential idea is to replace the traditional discrete representation of CPO by a more economical continuous representation in terms of a small number of analytical functions that we call structured basis functions (SBFs.) The number of SBFs required is equal to the number S of active slip systems in the crystals composing the aggregate (S = 3 for olivine), and each SBF represents the virtual CPO that would be produced by the action of just one of those systems. A CPO produced by the conjoint action of all slip systems is then represented as a superposition of the SBFs, with deformation-dependent coefficients chosen to give the best fit to the predictions of a homogenization model such as SOSC.

In addition to a continuous representation of CPO, our model involves a second important simplification: we assume that the CPO is a unique function of the finite strain experienced by the aggregate. Since many different deformation pathways can produce the same finite strain ellipsoid (FSE), our assumption implies that the CPO is independent of the deformation path (supposing of course that the FSE is initially spherical, corresponding to zero strain). The hypothesis of a unique relation between CPO and finite strain was first proposed by [START_REF] Mckenzie | Finite deformation during fluid flow[END_REF]. [START_REF] Ribe | On the relation between seismic anisotropy and finite strain[END_REF] verified this hypothesis using 2-D analytical and 3-D numerical models in which the crystallographic spin was a unique function of the crystal's orientation and did not depend on the overall orientation distribution. He also proposed an analytical parametrization of the percent of V p and V s anisotropy in the three principal planes of the FSE. [START_REF] Blackman | Seismic anisotropy in the upper mantle 2. Predictions for current plate boundary flow models[END_REF] compared VPSC and finite-strain predictions of CPO for 2-D and 2.5-D flow models of ridges and subduction zones, and found that the two methods yielded acceptably similar results over most of the flow field. More recently, [START_REF] Li | Subduction-induced mantle flow, finite strain, and seismic anisotropy: Numerical modeling[END_REF]) tested McKenzie's (1979) finite strain hypothesis using the VPSC model to calculate CPO along complex deformation paths in unsteady 3-D flows surrounding subducting slabs. Their figs 4-7 show that the peak orientations of the CPO are in general well predicted by the orientations of the axes of the FSE. For the particular values of the relative strengths of the three dominant slip systems assumed by [START_REF] Li | Subduction-induced mantle flow, finite strain, and seismic anisotropy: Numerical modeling[END_REF], the crystallographic axes [100], [010] and [001] are aligned with the long, short and intermediate axes of the FSE, respectively. Numerous authors have used the finite strain hypothesis to predict and/or interpret measurements of seismic anisotropy; examples include Phipps Morgan (1987), [START_REF] Silver | Shear wave splitting and subcontinental mantle deformation[END_REF], [START_REF] Hall | The influence of plate motions on three-dimensional back arc mantle flow and shear wave splitting[END_REF], [START_REF] Becker | Comparison of azimuthal seismic anisotropy from surface waves and finite strain from global mantle-circulation models[END_REF] and [START_REF] Long | Twodimensional modelling of subduction zone anisotropy with application to southwestern Japan[END_REF].

Our analytical representation of CPO together with the finite strain hypothesis yield an algorithm that is faster than standard kinematic models by a factor of ≈10 3 , and faster than homogenization approaches by a factor of ≈10 7 . It is also very simple to implement, as all the difficult calculations are done in advance. All that remains for the modeller to do is to follow the accumulation of finite strain along particle paths. In short, the simplicity and extreme rapidity of the present model makes it suitable for the calculation of CPO along pathlines in time-dependent flows of any complexity.

Our study is organized as follows. In Section 2, we review the Eulerian angles (φ, θ, ψ) ≡ g that describe the orientation of individual crystals; the orientation distribution function (ODF) f (g, t) that describes distributions of those orientations; and the evolution equation satisfied by f (g, t). In Section 3, we use a simple minimization principle to derive an analytical expression for the orientation-dependent crystallographic spin ġ(g) of a crystal with a single slip system in response to an externally imposed triaxial straining deformation. The resulting expression for ġ(g) is validated in Section 4, where it is shown to agree exactly (to within a multiplicative constant) with the predictions of the SOSC model. In Section 5, we solve the evolution equation for f (g, t) together with the analytical expression for ġ(g) to obtain the SBFs for the three dominant slip systems of olivine. Section 6 derives a quantitative parametrization (model SBF-1) for the coefficients in the SBF expansion of the CPO as functions of the axial ratios of the FSE and the relative strengths of the three slip systems. Section 6.2 presents an alternate parametrization (model SBF-2) in which ġ(g) is modified so as to predict textures that are in better agreement with experimental data than those predicted by the SOSC model. Finally, Section 8 proposes a simple modification of the calculation of the FSE that prevents singular textures from developing at large strain, and summarizes the advantages and limitations of our new parametrization.

D E S C R I B I N G C RY S TA L O R I E N TAT I O N A N D O R I E N TAT I O N D I S T R I B U T I O N

Our object of study is a polycrystalline aggregate comprising some large number N of olivine crystals. When the aggregate as a whole is subject to a given macroscopic deformation, its constituent crystals respond by deforming via internal shear on a small number S of slip systems. Each slip system s = 1, 2, ..., S is defined by a unit vector n s normal to the slip plane and a unit vector l s parallel to the slip direction.

The orientation of each crystal can be described by specifying three Eulerian angles that describe the crystal orientation relative to fixed external axes. In the texture community the most commonly used convention is that of [START_REF] Bunge | Texture Analysis in Materials Science[END_REF], which refers to consecutive rotations about the [001], [100] and [001] axes by angles (ϕ 1 , , ϕ 2 ). These angles are shown in Fig. 1 (a). In this study, however, it proves convenient to use a different set of Euler angles that are adapted to a slip system with slip direction l and unit normal n to the slip plane (Fig. 1b). In this convention the consecutive rotations are carried out around the axes l × n, l and l × n by angles (φ, θ, ψ). The orientation of the axes l, n and l × n relative to a sample coordinate system can also be characterized by the matrix of direction cosines

a i j (g) = ⎛ ⎜ ⎝ cφ cψ -sφ sψ cθ sφ cψ + cφ sψ cθ sψ sθ -cφ sψ -sφ cψ cθ -sφ sψ + cφ cψ cθ cψ sθ sφ sθ -cφ sθ cθ ⎞ ⎟ ⎠ , (1)
where 'c' and 's' indicate the cosine and sine, respectively, of the angle immediately following. Note that this new Euler angle convention is defined such that both n and l lie in the plane perpendicular to the axis that corresponds to changes in the angle ψ. Physically, this means that the crystallographic spin ġ produced by slip on the system has only a single non-zero component ψ; the other two components φ and θ are identically zero. This results in a significant simplification of the equation governing the evolution of the crystal orientation distribution, as we shall see below.

In the limit as the number of grains N → ∞, the distribution of their orientations can be described by a continuous function f, the orientation distribution function (ODF). In the case of texture evolution we consider the ODF as a time-dependent function f (g, t) defined such that f (g, t)dg is the volume fraction of crystals with orientations in an element dg of the Euler space centred on g at time t.

For crystals with the lowest possible (triclinic) symmetry, the ranges of the Eulierian angles required to describe a crystal's orientation are ϕ 1 ∈ [0, 2π ], ∈ [0, π ] and ϕ 2 ∈ [0, 2π ]. However, for crystals with orthorhombic symmetry such as olivine, the ranges are reduced to ϕ 1 ∈ [0, π ], ∈ [0, π ] and ϕ 2 ∈ [0, π ]. Using a proper normalization, the ODF representing a uniform (random) distribution is f = 1. Accordingly, the anisotropic part of an arbitrary ODF f is

f = f -1.
(2)

Throughout this paper, all ODFs and associated pole figures are expressed as multiples of a random distribution.

The time evolution of the ODF is governed by the continuity equation [START_REF] Clement | Prediction of deformation texture using a physical principle of conservation[END_REF])

0 = ∂ f ∂t + ∇•( ġ f ) ≡ ∂ f ∂t + ∂ ∂φ ( φ f ) + 1 sin θ ∂ ∂θ ( θ sin θ f ) + ∂ ∂ψ ( ψ f ), (3)
where ( φ, ψ, θ ) ≡ ġ is the rate of change of the orientation of an individual crystal. Eq. ( 3) is a conservation law that states that the rate of change of the volume fraction of crystals having orientations in a small element dg of the Euler space is equal to the net flux of crystal orientations into this element.

In general, the rotation rate ġ of a crystal depends on the imposed macroscopic strain rate tensor E ij , the crystal's orientation g, the overall orientation distribution f and the time t.

S I N G L E -S L I P M O D E L

The crystallographic spin for isotropic ODFs

Our next step is to obtain an explicit analytical expression for the crystallographic spin ġ as a function of the crystal orientation g and the imposed macroscopic velocity gradient tensor. In general, the velocity gradient tensor is the sum of a symmetric strain rate tensor E ij and an antisymmetric rotation rate tensor ij that measures the vorticity of the flow. However, because we assume that texture depends only on finite strain, we can set ij = 0 without the loss of generality. The deformation then has the form of irrotational triaxial straining, and is completely described by the strain rate tensor E ij .

We consider the idealized case of an isotropic (f = 1) aggregate of crystals each of which has only a single-slip system. These idealizations seem at first sight to be highly restrictive and unrealistic, but this is not the case. It turns out that when all three slip systems of olivine are active together, as they are in the SOSC model, the orientation-dependence of the slip rate on each system considered separately is exactly that predicted by the single-slip model we are about to describe. Moreover, this exact agreement is found to persist as the texture evolves beyond its initial isotropic state. The strong assumptions of the single-slip model are therefore entirely appropriate.

To specify the model, let n(g) be a unit vector normal to the slip plane and l(g) a unit vector in the slip direction for a crystal with orientation g. When the aggregate is deformed, each crystal within it responds by deforming in simple shear on planes normal to n(g) at a rate γ (g). Accordingly, the local strain rate tensor inside the crystal is

e i j = γ S i j , ( 4
)
where

S i j = 1 2 (n i l j + n j l i ) (5)
is the Schmid tensor. Following [START_REF] Ribe | A theory for plastic deformation and textural evolution of olivine polycrystals[END_REF], we determine γ (g) by minimizing the integrated misfit

E = 1 2π 2 e i j -E i j e i j -E i j dg (6)
between the local and global strain rate tensors, where the integral is over the whole Euler space and repeated indices are summed. The factor 1/(2π 2 ) is added because

dg = 2π 2 . ( 7
)
The minimization of (6) can be performed either directly or subject to the constraint of global strain rate compatibility, viz.,

1 2π 2 e i j dg = E i j . ( 8
)
The result is (Appendix A)

γ = 2AS i j E i j , (9)
where A = 5 if the constraint (8) is imposed and A = 1 if it is not (i.e. if the minimization is done for each crystal separately.) The large difference between the two values of A is due to the fact that a single slip system can satisfy the global strain rate compatibility condition (8) by itself only if it slips very rapidly. Crystals with three slip systems have more degrees of freedom and consequently can satisfy global strain rate compatibility with smaller values of A. This is demonstrated in Appendix A for a 'triple-slip' generalization of our single-slip model, which predicts A = 5/3 if all three slip systems have the same strength. In the next section we confirm this result by predicting numerically that A = 1.663.

Returning to the single-slip model, we note that the crystallographic spin associated with the slip rate γ is

ω i (g) = 1 2 i jk n j (g)l k (g) γ (g), ( 10
)
where ε ijk is the Levi-Civita permutation tensor. The corresponding rates of change of the Eulerian angles are

φ = ω 3 + (ω 2 cos φ -ω 1 sin φ) cot θ, (11a) θ = ω 1 cos φ + ω 2 sin φ, (11b) ψ = (ω 1 sin φ -ω 2 cos φ) csc θ. ( 11c
)
The final step is to write the spins φ, θ and ψ defined by ( 10) and ( 11) explicitly in terms of the Eulerian angles (φ, ψ, θ). We suppose without loss of generality that the external reference frame coincides with the principal axes of E ij . Let the principal strain rates be E 11 ≡ ˙ 1 , E 22 = ˙ 2 and E 33 ≡ ˙ 3 , where ˙ 1 + ˙ 2 + ˙ 3 = 0 because the deformation is incompressible. Using the slip system-specific Eulerian angles (φ, ψ, θ) defined in Fig. 1(b), we see that slip in the direction l on planes normal to n produces a crystallographic spin in the l × n direction. Therefore φ = θ = 0 identically, leaving ψ as the only non-zero spin component. Fig. 1(b) implies that l i = a 1i and n i = a 2i , where a ij are the direction cosines (1). Eqs ( 10) and ( 11) then become

φ = θ = 0, ψ = A(P cos 2ψ + Q sin 2ψ), (12a)
P = 1 2 (˙ 1 -˙ 2 )F(φ, θ), (12b) Q = 1 2 [(˙ 1 -˙ 2 )G(φ, θ) + (˙ 2 -˙ 3 )H (θ )] , (12c) F(φ, θ) = -sin 2φ cos θ, (12d) G(φ, θ) = sin 2 φ cos 2 θ -cos 2 φ, ( 12e
)
H (θ) = -sin 2 θ. ( 12f
)
In closing this subsection, it may be helpful to compare our single-slip model with the well-known models for polycrystal plasticity of [START_REF] Sachs | Zur Ableitung einer Fliessbedingung[END_REF] and [START_REF] Taylor | Plastic strain in metals[END_REF]. The Sachs model assumes that the stress within all grains is the same, and that only a single-slip system (the one with the most favourable Schmid factor) is active. While intergrain stress equilibrium is ensured, strain rate compatibility is not, with the consequence that voids and overlaps arise between the grains. By contrast, the Taylor model assumes that all grains experience the same strain rate, and that slip occurs on multiple slip systems simultaneously. The Taylor model is in a sense the opposite of the Sachs model: intergrain strain rate compatibility is ensured, but stress equilibrium is violated. Our model is similar to the Sachs model in that only one slip system is active, but differs from it in that intergrain stress equilibrium is not satisfied -in fact, the stress within the grains is not even calculated. Our model is similar to Taylor's in that strain rate compatibility is satisfied in a global sense [assuming that (8) is enforced], but differs from it in that different grains experience different strain rates depending on their orientations.

Structured basis functions: solving the evolution equation

Our goal in this subsection is to obtain an analytical solution of the evolution eq. ( 3) with the expressions (12) for the spins φ, θ and ψ. As in the previous section we consider the case of triaxial straining with strain rates ˙ 1 , ˙ 2 and ˙ 3 . Supposing that these strain rates are constants, the logarithms of the ratios of the axes c 1 , c 2 and c 3 of the FSE are

r 12 ≡ ln c 1 c 2 = (˙ 1 -˙ 2 )t, r 23 ≡ ln c 2 c 3 = (˙ 2 -˙ 3 )t, r 31 ≡ ln c 3 c 1 = -r 12 -r 23 .
For later convenience we also define the equivalent strain r 0 as

r 0 = 2 3 (˙ 2 1 + ˙ 2 2 + ˙ 2 3 ) 1/2 t = 2 3 r 2 12 + r 12 r 23 + r 2 23 1/2 . ( 13
)
Although we derived (12) under the assumption of an isotropic ODF we shall see in Section 4 that it is valid also for the anisotropic case if we allow the amplitude A to depend on the strains r 12 and r 23 . Since strain is a function of time, so also is A.

For solving the evolution equation it will prove useful to introduce an enhanced time

T (t) = t 0 A((˙ 1 -˙ 2 )τ, (˙ 2 -˙ 3 )τ )dτ. ( 14
)
Since ∂/∂t = A∂/∂T, the amplitude A cancels out of the evolution equation when the latter is written in terms of T. Based on the enhanced time we also define enhanced finite strain parameters

R 12 = (˙ 1 -˙ 2 )T (t), R 23 = (˙ 2 -˙ 3 )T (t), R 31 = (˙ 3 -˙ 1 )T (t), (15)
or equivalently,

˙ 1 T = 1 3 (R 12 -R 31 ), ˙ 2 T = 1 3 (R 23 -R 12 ), ˙ 3 T = 1 3 (R 31 -R 23 ). ( 16
)
When A = 1, T reduces to the true time t, while R ij reduce to the true (logarithms of the) axial ratios r ij of the FSE. Note that R 12 + R 23 + R 31 = 0 identically because the deformation is incompressible.

Because the only non-zero component of the crystallographic spin is ψ, the evolution eq. ( 3) reduces to

∂ f ∂t + ∂ ∂ψ ( ψ f ) = 0, ( 17
)
where ψ is given by ( 12). The above equation can be solved by the method of characteristics subject to the initial condition

f(t = 0) = 1 (Appendix B). The result is f (t, φ, θ, ψ) = 1 cosh( P2 + Q2 ) 1/2 -sin 2(ψ -χ ) sinh( P2 + Q2 ) 1/2 ≡ F, ( 18
)
where

P = 1 2 R 12 F(φ, θ), Q = 1 2 [R 12 G(φ, θ) + R 23 H (θ )] ( 19
)
and tan 2χ = Q/ P. In the rest of this paper we denote the analytical solution (18) by the special symbol F.

The solution F plays a fundamental role in this study: it is a 'structured basis function' that serves as a building block to represent complex textures that result from the simultaneous activity of several slip systems. The qualifier 'structured' indicates that the SBF automatically satisfies the symmetries associated with the slip system geometry and the imposed deformation, unlike typical unstructured basis functions such as generalized spherical harmonics. This means that the SBF conforms to the shape of an evolving texture. Moreover, the SBF automatically becomes more concentrated as the deformation progresses, as real textures do. Note that the SBF depends on the enhanced strains R 12 and R 23 rather than on the true strains r 12 and r 23 .

Because the general SBF ( 18) is a function of the slip system-specific Eulerian angles, there are in fact three distinct SBFs, one for each dominant slip system. These three SBFs, denoted by the symbols F s (s = 1, 2, 3), are related to one another by simple rotation matrices.

Let n and l be the normal and the direction of a slip system and f(t, φ, θ, ψ) the solution (18) of the evolution equation with respect to the slip system-related Euler reference frame. In order to evaluate f at an orientation g given by Bunge Euler angles (ϕ 1 , , ϕ 2 ) we proceed as follows:

(1) Determine a rotation M such that M z = n × l and M x = l.

(2) Determine the rotation matrix R corresponding to the Bunge Euler angles (ϕ 1 , , ϕ 2 ).

(3) Determine the Bunge Euler angles ( φ1 , ˜ , φ2 ) corresponding to the matrix M R. (4) The Bunge Euler angles ( φ1 , ˜ , φ2 ) just coincide with the slip system-based Euler angles (φ, θ, ψ) and can be directly plugged into (18).

As a consequence, the SBF with respect to any arbitrary slip system is just a rotated version of the ODF f defined in (18), namely, by the matrix M.

We now turn to the specific case of olivine with its three dominant slip systems (010)[100], (001)[100] and (010) [001]. Let the corresponding SBFs given by (18) be F 1 , F 2 and F 3 , respectively. These basis functions depend only on the orientation g and the slip system-specific enhanced strains

R (s) i j = (˙ i -˙ j )T s , T s (t) = t 0 A s (τ )dτ, ( 20
)
where A s (s = 1, 2, 3) are the (in general different) spin amplitudes for the three slip systems. As an illustration, Fig. 2 shows the pole figures of F 1 , F 2 and F 3 for uniaxial compression along the z-axis with enhanced strains R As a second illustration, Fig. 3 shows the pole figures of the SBFs F 1 , F 2 and F 3 for plane strain with enhanced strains R

(s) 12 = R (s) 23 = 1 and R (s) 31 = -2.
The pole figures are no longer symmetric about the z-axis due to the lower symmetry of plane strain compared with uniaxial compression.

As we noted above, the SBFs F s become progressively sharper as the strain increases. We illustrate this for F 1 and a deformation in uniaxial compression along the z-axis with enhanced finite strain parameters R A final important point is that the functions F 1 , F 2 and F 3 are linearly dependent to lowest order in the small-strain limit R ij → 0. The forms of Fs = F s -1 in this limit can be obtained from ( 18), but a more revealing approach is to solve directly the evolution equation

∂ t F s + ∇• ( ġs F s ) = ∂ t Fs + ∇• ġs 1 + Fs = 0, ( 21
)
where the repeated index s is not summed. In the limit R ij → 0, 1 + Fs ≈ 1, whence ( 21 A direct calculation using the explicit expressions for ġs (Appendix C) shows that F1

A 1 (0) - F2 A 2 (0) - F3 A 3 (0) = 0 for |R (s) i j | 1, ( 23
)
where A s (0) are the (generally different) spin amplitudes for the slip systems s = 1, 2 and 3 at the initial time t = 0. Eq. ( 23) implies that F 1 , F 2 and F 3 are linearly dependent in the small-strain limit |R (s)

i j | 1. Now A 1 (0) = A 2 (0) = A 3 (0)
for the special case of slip systems with equal strength, whence (23) becomes lim

R (s) i j →0 (F 1 -F 2 -F 3 ) = -1.
(24) Fig. 5 shows the pole figures of

F 1 -F 2 -F 3 for uniaxial compression with A 1 (0) = A 2 (0) = A 3 (0), R (s) 12 = 0, and R (s)
= 0.15. The [100] and [010] pole figures are exactly equal to -1 everywhere, in accord with (24). The [001] pole figure departs from -1 by a small amount (±2 per cent) because the enhanced strain R (s) 23 = 0.15 is not small enough for (24) to hold exactly (we chose it that way for purposes of illustration).

In closing, we note that the only unknown parameters in the SBFs are the spin amplitudes A s (r 12 , r 23 ) for the three slip systems, which appear in the definition (14) of the enhanced time. We now determine these amplitudes by comparing (12) with the predictions of the SOSC model.

C O M PA R I S O N W I T H T H E S E C O N D -O R D E R S E L F -C O N S I S T E N T M O D E L

The single-slip model is based on a minimization principle that is purely ad hoc, and which therefore requires more rigorous justification. To this end, we now compare the predictions of the single-slip model with those of a more complete polycrystal plasticity model, the SOSC model of Ponte [START_REF] Ponte Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations: I-theory[END_REF]. The SOSC model is the most sophisticated member of the class of homogenization models discussed briefly in the Introduction. In most versions of the homogenization scheme, including the VPSC model, each mechanical phase is equivalent to a single grain, and the local (grain-scale) stress and strain rate tensors vary among the grains but not within them. The SOSC model goes a step further by treating a phase as an ensemble of grains with the same orientation and material properties, and allowing the stress and strain rate to vary both among and within the phases. As a result, the SOSC model predictions of quantities such as the effective average stress in a deforming aggregate agree significantly better with full-field calculations than those of competing homogenization schemes (such as VPSC) that do not account for intraphase stress and strain heterogeneity. The treatment of rheological nonlinearity in the SOSC model is also superior to that in the VPSC model, which suffers from a theoretical inconsistency first pointed out by [START_REF] Masson | An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals[END_REF].

The free parameter in our single-slip model is the amplitude A, which we aim to define in this section as a function of CRSS ratios p 12 and p 23 (defined hereafter) and the axial ratios r 12 and r 23 of the FSE. By calibrating A to fit the predictions of the SOSC model, we are able to achieve very high accuracy (typically > 99.5 per cent variance reduction). It is surprising that the amplitude A can be chosen independently of the grain orientation and of the global ODF and still provide an excellent fit to the SOSC predictions. This serves as a validation of our model, which is much simpler than the SOSC model and involves far fewer parameters.

The SOSC model requires that each crystal in the aggregate satisfy von Mises's criterion, according to which a crystal can only accommodate an arbitrary imposed deformation if it has at least five independent slip systems (although four independent slip systems are often found to be adequate in numerical practice; see [START_REF] Detrez | Effective viscoplastic behavior of polycrystalline aggregates lacking four independent slip systems inferred from homogenization methods; application to olivine[END_REF]. We assume that each olivine crystal has the three dominant slip systems (010)[100] (s = 1), (001)[100] (s = 2) and ( 010)[001] (s = 3), plus two harder slip systems (101)[-1 0 1] (s = 4) and (1 0 -1)[1 0 1] (s = 5). In the SOSC model, the slip rate γs on each slip system s obeys a power-law rheology of the form

γs ∝ τ τ s ms -1 τ τ s , ( 25
)
where τ is the resolved shear stress (i.e. the shear stress acting on the slip plane in the slip direction), τ s is a critical resolved shear stress (CRSS) that measures the inherent resistance of the slip system to slip, and m s is a power-law exponent. We assume m s = 3.5 for all slip systems [START_REF] Bai | High-temperature creep of olivine single crystals, 1. Mechanical results for buffered samples[END_REF]. Because the macroscopic deformation rate of the aggregate is specified in our SOSC calculations, only the ratios of the parameters τ s (and not their absolute values) are relevant. In our calculations we assume τ 1 /τ 2 ∈ [0.25, 4.0], τ 2 /τ 3 ∈ [0.25, 4.0] and τ 4 /τ 2 = τ 5 /τ 2 = 100 (Table 1). In the sequel, we shall characterize the CRSS ratios of the dominant slip systems s = 1, 2 and 3 in terms of the logarithmic ratios

p 12 = ln τ 1 τ 2 , p 23 = ln τ 2 τ 3 . ( 26
)
The comparison of the SOSC and single-slip model predictions proceeds as follows. First, the type of macroscopic deformation experienced by the aggregate (uniaxial compression, plane strain, etc.) is chosen by specifying ˙ 1 , ˙ 2 and ˙ 3 . Next, the initial condition for the SOSC model is specified; we use an ensemble of 2000 olivine orientations with an isotropic (random) distribution. The SOSC code is then run to determine the evolving orientations g (k) (t) and slip rates γs (g (k) , t) for each phase k and each slip system s. The calculated values of γs (g (k) (t)) are then substituted into (10) to obtain the partial spins ω (s) i (g (k) (t)) due to the actions of the individual slip systems, which are related to the total spin ω i by

ω i = 3 s=1 ω (s) i . ( 27
)
In writing ( 27), we have ignored slip systems s = 4 and 5 that make negligible contributions to the rotation. Finally, by substituting the partial spins ω (s)

i into (11) and expressing the results in terms of the Eulerian angles, we obtain the rotation rates ψs (g (k) s ) for all orientations k and slip systems s, which can be compared directly with the prediction (12) of the single-slip model.

Since the SOSC model involves three dominant slip systems, the predicted slip rates γs for each of those systems can be compared individually with the single-slip model. Let A 1 , A 2 and A 3 be the unknown amplitudes of the slip rates for each of the three dominant slip systems. As an example, we will show how the amplitude A 1 is determined. Fig. 6 shows the SOSC predictions of γ1 as a function of θ and ψ for two degrees of finite deformation in uniaxial compression:

|˙ 3 |t = 0 [part (a)] and |˙ 3 |t = 0.4 [part(c)].
To compare these functions with the predictions of the single-slip model, we evaluate the expression ( 9) with E 11 = E 22 = ˙ 0 /2 and E 33 = -˙ 0 , the values corresponding to uniaxial compression. Eq. ( 9) then reduces to

γ1 = 3 4 A 1 ˙ 0 sin 2 θ sin 2ψ. ( 28
)
Fitting (28) to the functions shown in Figs 6(a) and (c) by adjusting A 1 , we find A 1 = 1.66 and A 1 = 2.28, respectively. The corresponding functions (28) are shown in panels (b) and (d), and are indistinguishable from panels (a) and (c), respectively. A more quantitative comparison shows that the variance reduction (VR) is 99.9 per cent in both cases. This result, which was obtained by adjusting only a single parameter, validates the simple prediction of the single-slip model and shows that it holds for each slip system individually even when several slip systems are active together. An important corollary is that the angular dependence of γs remains the same as strain accumulates; only the amplitude A s of the function changes.

In the most general case, the amplitudes A s are functions of the CRSS ratios p 12 and p 23 and the axial ratios r 12 and r 23 of the FSE, that is, However, it turns out that to an excellent approximation the amplitudes A s are products of a function of p 12 and p 23 and a function of r 12 and r 23 , viz., A s = 1.663H s ( p 12 , p 23 )J (r 12 , r 23 ) where H s (0, 0) = J (0, 0) = 1.

A s = A s ( p 12 , p 23 , r 12 , r 23 ) (s = 1, 2, 3). ( 29
(

) 30
The functions H s describe how the slip system activities depend on the CRSS ratios at the initial instant (t = 0) of the deformation. The function J describes how the activities of slip systems with equal strengths vary as functions of finite strain for different types of deformation (uniaxial compression, plane strain, etc.). The prefactor 1.663 in (30) is the spin amplitude at the initial instant of the deformation (r 12 = r 23 = 0) when the three slip systems have equal strengths (p 12 = p 23 = 0). In fact, its value can be determined analytically using an extended 'triple slip' version of the minimization principle of Section 3 that allows all three slip systems of olivine to act simultaneously (Appendix A). The resulting value of the prefactor is exactly 5/3, which agrees to within 0.2 per cent with the value 1.663 determined numerically by fitting the slip rate (9) to the SOSC model predictions.

Our task now is to determine the functions H s ( p 12 , p 23 ) and J (r 12 , r 23 ) that appear in (30). Consider first the function J , which corresponds to p 12 = p 23 = 0. As illustrations, we consider two simple deformations: uniaxial compression (˙ 3 = -˙ 0 , ˙ 1 = ˙ 2 = ˙ 0 /2), and plane strain (˙ 1 = -˙ 3 = ˙ 0 , ˙ 2 = 0). Fig. 7 shows the best-fitting amplitudes A 1 -A 3 for these two types of deformation, as functions of the equivalent strain r 0 . Two features of these results are noteworthy. First, A 1 = A 2 = A 3 to within an error that ranges from 0.1 per cent for small strains to a few percent for large strains. Second, the VR of the fit exceeds 99 per cent in all cases for r 0 < 0.75. This implies that the single-slip model gives a nearly perfect fit to the SOSC model when the amplitudes A 1 -A 3 are chosen correctly.

We next determine the amplitudes A 1 = A 2 = A 3 for arbitrary finite strain parameters r 12 and r 23 , still assuming p 12 = p 23 = 0. As a representative value of the amplitude, we use the average of those values among A 1 , A 2 and A 3 for which the VR of the fit to the SOSC model is at least 99.5 per cent. It is convenient to plot the resulting amplitude on a so-called Flinn diagram with r 23 as the horizontal coordinate and r 12 as the vertical coordinate. On such a diagram the horizontal axis r 12 = 0 corresponds to FSEs produced by uniaxial compression, the vertical axis r 23 = 0 to FSEs produced by uniaxial extension and the line r 12 = r 23 to FSEs produced by plane strain (pure shear). Fig. 8(a) shows the amplitude of the single-slip model that best fits the predictions of the SOSC model. To make the connection with Fig. 7, note that the curves for uniaxial compression and plane strain in that figure are cuts along the horizontal axis r 12 = 0 and the diagonal r 12 = r 23 , respectively, of the contour plot of Fig. 8(a). To use the results of Fig. 8(a) in a numerical code, we need to fit them using an analytical function J . The explicit form of this function is given in Appendix D, and a contour plot of the quantity 1.663J is shown in Fig. 8(b).

Turning next to the limit r 12 = r 23 = 0, we again use a least-squares approach to determine the amplitudes A 1 -A 3 of the single-slip rotation rate (12) that best fit the SOSC predictions. In all cases the fit is nearly perfect, with a variance reduction VR ≥ 99.7 per cent. The last amplitude-related quantity that we shall need is the enhanced time T s defined by ( 14). More precisely, we are interested in the quantity T s /t, which is just the ratio of the enhanced finite strain parameters R where H s is the same function that appears in (30). The function K is determined by numerical integration of the amplitude shown in Fig. 8(b), and the result is shown in Fig. 10(a). An explicit analytical expression for K(r 12 , r 23 ) is given in Appendix D, and a contour plot of it is shown in Fig. 10(b).

To summarize: by fitting the expression (12) for the crystallographic spin to the predictions of the SOSC model, we have obtained analytical expressions for the spin amplitudes A s and the related quantity t -1 T s . We call this set of analytical expressions 'spin model 1'. It is defined by

A s = 1.663H s ( p)J (r), t -1 T s = H s ( p)K(r) (spin model 1), (32)
where p = ( p 12 , p 23 ), r = (r 12 , r 23 ) and the explicit expressions for H s , J and K are given in Appendix D. In spin model 1, the amplitudes of the crystallographic spins increase strongly with increasing r. However, we shall also use a simpler 'spin model 2' defined by

A s = t -1 T s = H s ( p), (spin model 2). ( 33
)
In spin model 2, the amplitudes of the spins are independent of the strain r. The motivation for this second model is explained below in Section 6.2.

N U M E R I C A L S O L U T I O N O F T H E E V O L U T I O N E Q UAT I O N

Now that we have explicit analytical expressions for the spin amplitudes A s , we no longer need to use the computationally intensive SOSC model. Instead, we can proceed much more efficiently by solving numerically the evolution eq. ( 3). There are two main advantages to this strategy. First, the SOSC model does not deliver the ODF as such, but only an ensemble of orientations that must be transformed into an ODF by means of certain mathematical approximations. By contrast, numerical solution of (3) yields the ODF directly. The second advantage is that we can obtain solutions for the ODF using expressions for the spin amplitudes A s that differ from those predicted by the SOSC model. One of these is the spin model 2 that was briefly mentioned at the end of the previous section.

The evolution eq. ( 3) can be solved numerically by discretization in either Fourier space or in orientation space. Discretization in Fourier space works especially well when the deformation is small, resulting in smooth textures that can be well approximated by generalized spherical harmonics of low degree [START_REF] Barton | The use of discrete harmonics in direct multi-scale embedding of polycrystal plasticity[END_REF]. However, we seek a model that can handle large deformations and sharp textures, and so we choose to discretize the orientation space. We do so by specifying an initially uniformly spaced grid of 5832 orientations ('grains'), and we employ the method of characteristics to solve the evolution eq. ( 3) subject to an isotropic initial condition f = 1 for all orientations.

As deformation proceeds, the orientation of each grain will wander in the Euler space at a rate and in a direction governed by the inherent crystallographic spin for each orientation g. The paths along which the grain orientations wander are called the characteristic curves. Each curve is governed by the first-order differential equations

dφ dt = φ(g), dθ dt = θ (g), dψ dt = ψ(g), (34)
where ( φ, θ, ψ) are the crystallographic spins with respect to the slip system-specific Eulerian angles (φ, θ, ψ). These spins are the sums of contributions from each of the three dominant slip systems of olivine, and are given in Appendix C. The amplitudes A s that appear in those expressions can be those of either spin model 1 or spin model 2. The eq. ( 34) is solved using a second-order Runge-Kutta (midpoint) method to determine a new point g(dt) on the characteristic curve that results from crystallographic spin operating during a small time interval dt. We then return to the original point (time t = 0) and use the midpoint method to evolve the ODF along the characteristic by solving the equation

d f dt = -f ∇ • 3 s=1 ġs (g) . ( 35
)
The process is then repeated for each characteristic curve, stepping the ODF forward by increments of dt until the desired maximum strain is reached. At this point, the grain orientations are no longer uniformly distributed in the Euler space: they are more concentrated where f is large and sparser where f is small.

E X PA N S I O N O F O R I E N TAT I O N D I S T R I B U T I O N F U N C T I O N S I N T E R M S O F S T RU C T U R E D B A S I S F U N C T I O N S

In the previous section we saw how to solve the evolution eq. ( 3) numerically for arbitrary CRSS ratios p and arbitrary strains r. We now show how the resulting ODF f num can be represented as a weighted sum of the SBFs for the three dominant slip systems of olivine. We posit the expansion

f num ≈ 1 + 3 s=1 C s ( p, r) Fs (g, R (s) ), ( 36
)
where

R (s) = (R (s) 12 , R (s)
23 ) and C 1 -C 3 are unknown coefficients that depend continuously on p and r. Because Fs are the anisotropic parts of the SBFs, the normalized integral of (36) over the Euler space is unity. Moreover, (36) automatically satisfies the isotropic initial condition f = 1 at r = 0.

For given values of p and r, the numerical solution of the evolution equation delivers an ODF f num (g k ) at N orientations g k (k = 1, 2, ..., N) in the Euler space. The coefficients C 1 -C 3 are then determined by minimizing the mean squared error where the arguments p, r and R (s) have been suppressed to simplify the notation. Setting to zero the derivatives of ( 37) with respect to C s yields three simultaneous linear equations for those coefficients that are solved numerically. However, this procedure cannot be used at zero strain (r = 0) because the functions Fs are not linearly independent in that limit (see 23). Accordingly, we determined the coefficients C s at zero strain by extrapolating backwards the values obtained for small equivalent strains r 0 = 0.01-0.02.

N k=1 f num (g k ) -1 - 3 s=1 C s Fs (g k ) 2 , ( 37
)
In the next two subsections, we use the procedure described above to determine the coefficients C 1 -C 3 for spin models 1 and 2.

Spin model 1

In this model, to repeat, the spin amplitudes A s = 1.663J (r)H s ( p) are those that match the predictions of the SOSC model. Fig. 11 shows the expansion coefficients C 1 -C 3 as functions of the equivalent strain r 0 for an aggregate with p 12 = p 23 = 0 deformed in either (a) uniaxial compression (˙ 1 = ˙ 2 = -˙ 3 /2) or (b) plane strain (˙ 2 = -˙ 1 , ˙ 3 = 0). The corresponding VR of the fit is shown by the dashed lines, according to the scale along the east edge of the figure. The coefficients for uniaxial compression and plane strain are identical at r 0 = 0 because the SBF expansion (36) satisfies the evolution eq. (3) exactly in the limit r 0 → 0. The VR of the fit is therefore 100 per cent for r 0 = 0. The VR decreases as r 0 increases, because the expansion ( 36) is no longer an exact solution of the evolution equation. For r 0 ≤ 0.6, VR≥ 90 per cent for uniaxial compression and VR≥ 82 per cent for plane strain. Fig. 12 compares pole figures of the ODF obtained by the method of characteristics (top) with its best-fitting SBF expansion (bottom) for deformation in uniaxial compression to an equivalent strain r 0 = 0.27. The CRSS ratios are p 12 = p 23 = 0 in both cases. The two sets of pole figures are remarkably similar in overall pattern. However, there are minor differences in the amplitudes: the pole figures of the SBF expansion have somewhat lower maxima and minima.

Fig. 13 is analogous to Fig. 12, but for deformation in plane strain to an equivalent strain r 0 = 0.20. The systematics are the same as in Fig. 12: an overall similarity in pattern, with somewhat lower minima and maxima in the pole figures of the SBF expansion.

This is an appropriate place to remark that nothing about the assumed form of the expansion (36) prevents the ODF from taking on unphysical negative values over some portion of the Euler space. Indeed, explicit calculation shows that while (36) is always positive definite for small to moderate strains, it begins to exhibit negative values beyond some critical strain. The critical strains for Fig. 11 are r 0 = 0.275 for uniaxial compression and r 0 = 0.207 for plane strain, and are shown by small vertical black bars in each panel of the figure. However, ODFs with negative values are easy to correct by setting the negative values to zero and then rescaling the ODF so that its integral over the Euler space is unity. Dealing with spurious negative values of the ODF is a standard concern in texture analysis [START_REF] Bunge | Texture Analysis in Materials Science[END_REF].

Spin model 2

As we noted earlier, in spin model 2 the spin amplitudes A s = H s ( p) are independent of the strain r. To understand the motivation for this choice, it is necessary to examine pole figures for the texture obtained in a laboratory deformation experiment. Fig. 14(a) shows pole figures for an olivine aggregate experimentally shortened by 58 per cent in uniaxial compression (Nicolas et al. 1973). For comparison, Fig. 14(b) shows pole figures predicted by the SOSC model for 50 per cent shortening. Part (b) is for 50 per cent shortening rather than 58 per cent because the SOSC algorithm fails to converge beyond about 52 per cent shortening. The two sets of pole figures have the same general structure, but the ones predicted by the SOSC model are much more concentrated despite the smaller deformation. This indicates that the SOSC model predicts textures that are much too sharp to match laboratory deformation textures. One reason for this is that in reality, significant (≈10 • ) misorientations develop inside grains at large strain, an effect that tends to 'spread' the ODF and prevent very large values. As this orientation spread is not accounted for in the SOSC model, the latter predicts ODFs that are too sharp at strains typically larger than r 0 ≈ 0.5. However, the spreading effect by itself cannot account for the discrepancy between Figs 14(a) and (b), because Fig. 14(a) represents only non-recrystallized grains. Another likely reason is the strong requirement (8) of strain rate compatibility in the SOSC model, which forces the grains to spin ever faster as the deformation progresses. This acceleration is evident in Fig. 7, where it appears that the spin amplitudes A 1 -A 3 may become infinite at a finite value of the equivalent strain (finite-time singularity). In reality, strain rate compatibility is partly ensured by deformation in the recrystallized grains, and so we should not expect pole figures based solely on unrecrystallized grains to match the SOSC predictions.

A simple way to get around the problem just discussed is to refer to the derivation in Section 3, where we determined the spin to within an unknown amplitude A by minimizing the misfit between the macroscopic and microscopic strain rates. There we saw that the amplitude A

C s C s VR VR C 1 C 3 C 2 equivalent strain r 0 equivalent strain r 0 C 2 , C 3 C 1 Figure 15
. Same as Fig. 11, but for the model SBF-2.

= 1 (constant) if strain rate compatibility is not imposed as a constraint. Thus to obtain more reasonable textures, spin model 2 assumes A s = 1 for the case of slip systems having equal strengths. A simple generalization to the case of slip systems with different strengths is to set

A s = H s ( p).
The coefficients C 1 -C 3 for spin model 2 are determined in the same way as those for spin model 1. Fig. 15 shows the results for an aggregate with p 12 = p 23 = 0 deformed in either (a) uniaxial compression or (b) plane strain, and should be compared with the analogous Fig. 11 for slip model 1. Because the spin amplitudes A s do not increase with increasing strain in spin model 2, the evolution eq. ( 3) can be solved to large equivalent strains r 0 without difficulty, unlike the case of spin model 1 for which singular behaviour is encountered for moderate values of r 0 . The maximum value of r 0 in Fig. 15 is 0.9, as opposed to 0.6 in Fig. 11. A further advantage of spin model 2 is that the variance reduction VR of the SBF expansion is better than for spin model 1. In Fig. 15, for example, VR≥ 93.4 per cent for uniaxial compression and VR≥ 88.6 per cent for plane strain.

We also plotted pole figures analogous to Figs 12 and 13 with spin model 1 replaced by spin model 2 and the deformation carried out to larger equivalent strains r 0 = 0.47 and 0.34, respectively. However, we do not show those pole figures here because they are visually identical to Figs 12 and 13. The reason is that the values of r 0 in each figure are just slightly less than the critical values at which negative values start appearing in the SBF expansions. For each deformation type (uniaxial compression or plane strain), those critical values correspond to nearly the same values of the enhanced strains R (s) i j . Since the evolution of texture depends on R (s) i j rather than the true strains r ij , it is not surprising that the pole figures for spin models 1 and 2 display such strong similarities for appropriately chosen values of r 0 .

The use of spin model 2 rather than spin model 1 results in a much improved agreement with laboratory data. This is demonstrated by Fig. 14(c), which shows pole figures for 58 per cent shortening in uniaxial compression predicted using an SBF expansion with spin model 2 and the same values of p 12 and p 23 as for the SOSC model (Fig. 14b.) While the overall patterns of the SOSC and SBF pole figures are similar, the SBF amplitudes agree much better with the laboratory data (Fig. 14a) than do the SOSC amplitudes. For completeness, Fig. 14(d) shows the analogous pole figures predicted by the D-Rex model [START_REF] Kaminski | D-Rex, a program for calculation of seismic anisotropy in the convective upper mantle[END_REF]) without recrystallization. The maximum amplitude (≈6.5) is about 24 per cent less than that of the SBF model, and agrees somewhat better with the experimental pole figures of Fig. 14(a). However, the amplitudes of the pole figures of Figs 14(a),(b) and (d) should not be overinterpreted, because they depend on the value chosen for the half-width of the Gaussian function used to approximate the infinitely sharp (delta-function) orientation distribution of each individual grain. The smaller the half-width used, the sharper are the peaks in the pole figures. In Fig. 14(d), for example, a value = 7.5 • is used. If one instead chooses = 5 • , then the maximum amplitude of the pole figures increases to ≈7.5; and a still smaller value of would give a still larger amplitude. Given this fundamental indeterminacy of the amplitude, we conclude that the pole figures in Figs 14(c) and(d) are consistent with one another. Note that the maximum amplitude of the pole figures in Fig. 14(c) is unambiguous because those pole figures are based on fitting SBFs to a numerical solution (by the method of characteristics) of the evolution eq. ( 3). Because that solution delivers the ODF directly rather than a set of individual grain orientations, the choice of a value of is unnecessary.

Analytical expressions for the expansion coefficients

In general, the SBF expansion coefficients C s depend on the CRSS ratios p 12 and p 23 and on the strain parameters r 12 and r 23 . It turns out that the dependency can be described very accurately with a polynomial function of the form

C s ( p 12 , p 23 , r 12 , r 23 ) = 2 i=0 2-i j=0 4 k=0 4-k l=0 Ĉ(s) i jkl p i 12 p j 23 r k 12 r l 23 , s = 1, 2, 3, (38)
where i, j, k and l are powers. We determined the coefficients 386, 1.386], and r 0 ∈ [0, 0.6] (spin model 1) or r 0 ∈ [0, 0.9] (spin model 2) using a standard least-squares procedure. The final results are tabulated in Appendix E for spin model 1 and Appendix F for spin model 2. For spin model 1, the RMS errors of the fits to C 1 , C 2 and C 3 are 0.0223, 0.0135 and 0.0137, respectively. The corresponding errors for spin model 2 are 0.0164, 0.0108 and 0.0111.

Ĉ(s) i jkl for p 12 ∈ [ -1.386, 1.386], p 23 ∈ [ -1.

S U M M A RY A N D S T E P -B Y-S T E P G U I D E

For the reader's convenience we give here a brief summary of the method presented in this paper, followed by a step-by-step guide to its use.

Summary

(1) Currently available methods for calculating evolving CPO in upper mantle flow fields have a prohibitively high computational cost, which results in large part from an unwieldy discrete representation of CPO in terms of ≈10 3 -10 4 individual grains.

(2) As a first building block for a new approach, we present a 'single-slip' model for an aggregate of crystals each containing only a single-slip system and deformed by an externally imposed strain rate tensor E ij . By minimizing the misfit between E ij and the local (crystal) strain rate tensor, we predict that the rate of slip on slip planes is γ = AE i j S i j , where A is an amplitude factor and S ij is the Schmid tensor that depends on the orientation of the slip system. The derivation assumes that the ODF is isotropic.

(3) By comparing the above prediction of the slip rate with the prediction of the SOSC model with three active slip systems (s = 1, 2, 3), we find that the expression γ = A s E i j S i j is valid for each slip system individually even when the ODF is non-isotropic, as long as the amplitudes A s are allowed to be strain-dependent. The amplitudes A s are the same only when the CRSSs of the slip systems are the same.

(4) The result γ = AE i j S i j for the case of single slip allows the evolution eq. ( 3) for the ODF to be solved analytically for an aggregate of crystals having only a single slip system and being deformed by irrotational triaxial straining. The result is a 'structured basis function' that represents the virtual ODF that would be produced by a single slip system acting alone. There are three SBFs, one for each slip system, which are related to one another by finite rotations.

(5) We assume that CPO is a unique function of the finite strain experienced by the aggregate as long as the rate of change of the velocity gradient tensor along flow pathlines is not too fast. It is therefore sufficient to consider only deformation by triaxial straining, the kind of deformation assumed in the derivation of the SBFs.

(6) An ODF produced by the joint action of three slip systems can be represented to an excellent approximation by a weighted sum of the three SBFs. The coefficients C s (s = 1, 2, 3) in this expansion depend on the ratios p 12 and p 23 of the CRSSs of the slip systems, and on the ratios r 12 and r 23 of the axes of the FSE. We give analytical expressions for C s (p 12 , p 23 , r 12 , r 23 ) for two spin models: one in which A s increases with strain, and a second in which A s is strain-independent. The latter is our preferred model because it agrees better with the CPO obtained in laboratory experiments.

(7) Our new SBF expansion approach to CPO calculation is ≈10 3 times faster than kinematic models, and ≈10 7 times faster than homogenization models such as SOSC. Furthermore, a simple modification of the FSE enables the user to avoid unrealistic singular CPO at very large strains (see Concluding Remarks section).

Step-by-step guide

(1) Determine the orientation o and the axial ratios r 12 , r 23 of the FSE at some point of interest in the flow field.

(2) Choose the values of the critical resolved shear stress ratios p 12 and p 23 .

(3) Compute the SBF expansion coefficients (5) Calculate the three SBFs F s (g, R (s) ) where R (s) = (R

C O N C L U D I N G R E M A R K S

In this paper we use several different models to calculate CPO, so some further explanation of the differences among them may be useful.

The five models we use are summarized in Table 2 .

The most computationally intensive is the SOSC model which enforces compatibility of both stress and strain rate between each grain and a homogenous effective medium that represents the average of all the other grains. The model predicts amplitudes of the crystallographic spin that increase rapidly with strain, and is therefore limited to moderate equivalent strains r 0 ≈ 0.6. In the context of 3-D time-dependent convection models, the high computational cost of the SOSC model limits its use to a small number of selected pathlines. However, [START_REF] Blackman | Effects of crystal preferred orientation on upper-mantle flow near plate boundaries: rheologic feedbacks and seismic anisotropy[END_REF] recently showed that full coupling of the SOSC model to a steady 2-D flow field with texture-dependent anisotropic viscosity is possible.

A faster alternative to the SOSC model is either of the models CHAR-1 and CHAR-2, both of which use the method of characteristics to solve the evolution eq. ( 3) for the ODF. Model CHAR-1 uses analytical expressions for the spin amplitudes that match the predictions of the SOSC model. CHAR-2 uses spin amplitudes that do not increase with strain, and consequently predicts CPO in better agreement with experimental observations. For a given number of grains/characteristics, the CHAR-1 and CHAR-2 models run ≈10 4 times faster than the SOSC model. However, this is still too slow for routine use in geodynamics.

The (relatively) high computational cost of the CHAR-1 and CHAR-2 models motivates our final models, which we call SBF-1 and SBF-2. In these models, the ODF is expressed as a weighted sum of analytical SBFs, each of which represents the virtual CPO that would be produced by a single-slip system acting alone. Model SBF-1 (like CHAR-1) uses expressions for the spin amplitudes that match the predictions of the SOSC model (spin model 1), while model SBF-2 (like CHAR-2) assumes that the amplitudes do not increase with strain (spin model 2). The coefficients of the SBFs in the expansions are given analytically as functions of the CRSS ratios p 12 and p 23 and the strain parameters r 12 and r 23 . Both SBF-1 and SBF-2 are ≈10 3 times faster than models CHAR-1 and CHAR-2, and ≈10 7 times faster than SOSC. So that potential users can make a choice, we have provided the expansion coefficients for both the SBF-1 and the SBF-2 models in Appendices E and F, respectively. However, SBF-2 is our preferred model because it gives a better fit to experimental data.

In addition to their inherent computational speed, the CHAR and SBF models have another significant advantage over the SOSC model: they do not require the presence of fictitious slip systems. We saw that the SOSC model does not converge unless at least one additional hard slip system with a high CRSS is specified in addition to the three 'true' slip systems (010)[100], (001)[100] and (010) [001]. The CHAR and SBF models are not subject to this necessity, and allow one to focus entirely on the behaviour of the three true slip systems.

A further important advantage of our SBF models is that they provide a simple way to limit the sharpness of the predicted textures. It is well known that self-consistent models such as VPSC and SOSC predict textures that are too sharp relative to experimental observations, because the predicted amplitude of the crystallographic spin increases rapidly with increasing strain. Figs 14(b) and 16(b) showed examples of excessively sharp textures predicted by the SOSC model. To a lesser degree, the D-Rex model also predicts textures that are too sharp when the strain is very large, even though the spin amplitude in that model does not increase with strain. In our new SBF models, this undesirable effect can be avoided by adjusting the shape of the FSE so that it never exceeds a specified maximum equivalent strain r max 0 . Let r 0 be the true equivalent strain of the FSE at some point along a pathline, and let r 12 and r 23 be the corresponding logarithmic ratios of its axes. We propose to calculate the adjusted values of these parameters as

) 39
The formulae (39) correspond to moving the FSE along a radial line towards the origin in a Flinn diagram, a transformation that preserves the type of deformation (uniaxial compression, plane strain, etc.) encoded in the FSE. The formula for r adj 0 implies that r adj 0 = r 0 (i.e. no adjustment) for small strains and r adj 0 = r max 0 for large strains. The simple expedient just described means that in the limit of infinite (true) strain the texture reaches a steady state with a finite sharpness.

Our recommended model SBF-2 has a close relationship with the D-Rex model of [START_REF] Kaminski | D-Rex, a program for calculation of seismic anisotropy in the convective upper mantle[END_REF]. In fact, both models use exactly the same expression (12) for the crystallographic spin as a function of crystal orientation, with the same amplitudes A s = 1 when the slip systems have equal strengths. However, D-Rex is based on a discrete (many-grain) representation of the texture whereas the much faster SBF-2 model uses a continuous representation. D-Rex also predicts excessively sharp textures when the strain is large. On the other hand, D-Rex includes two factors that are missing in SBF-2, namely dynamic recrystallization and the presence of a second phase (enstatite). These factors are discussed in more detail below. It is also of interest to compare the models presented here with the ANPAR model of [START_REF] Goulding | Analytical parametrization of self-consistent polycrystal mechanics: Fast calculation of upper mantle anisotropy[END_REF]. The main difference is that ANPAR does not assume that texture is a unique function of finite strain. Accordingly, ANPAR allows the principal axes of the macroscopic strain rate tensor E ij to have an arbitrary orientation relative to the axes of the FSE, and provides an analytical expression for the instantaneous crystallographic spin for this very general case. This allows ANPAR to match with near exactitude the predictions of the SOSC model, even for rotational deformations such as simple shear. For the limiting case of irrotational deformation where the principal axes of E ij and the FSE coincide, ANPAR is equivalent to model CHAR-1 of this paper, and has a similar computational cost.

We conclude this section with a discussion of three limitations of our model, which are important to keep in mind. The first is the neglect of dynamic recrystallization. The influence of dynamic recrystallization on CPO development was demonstrated experimentally by [START_REF] Zhang | Lattice preferred orientation of olivine aggregates deformed in simple shear[END_REF], who deformed olivine aggregates in simple shear. In simple shear, the long axis of the FSE is initially at 45 • to the shear plane, and then rotates progressively towards that plane. [START_REF] Zhang | Lattice preferred orientation of olivine aggregates deformed in simple shear[END_REF] observed that the peak orientation of the fast (a-) axis of olivine followed the long axis of the FSE up to ≈ 50 per cent shear strain, but then rotated much faster towards the shear plane at larger strains. This is illustrated in Fig. 16(a), which shows pole figures for an initially isotropic olivine aggregate deformed in the laboratory in simple shear to a shear strain γ = 140 per cent [START_REF] Lee | A scanning electron microscope study of dynamic recrystallization on lattice preferred orientation in olivine[END_REF]. The orientation of the axis of maximum instantaneous extension rate is shown by the dashed line, and the orientation of the long axis of the FSE is shown by the solid line. The maxima of the [100] axis orientations are not in the direction of the long axis of the FSE, but instead lie close to the shear plane (red lines). For comparison, the remaining sets of pole figures in Fig. 16 show the predictions of three different models, all with p 12 = -0.693 and p 23 = -0.405. Fig. 16(b) shows the texture predicted by the SOSC model, which closely tracks the orientation of the FSE. Fig. 16(c) shows the pole figures predicted by the SBF-2 model. This texture is aligned with the FSE by definition, but has much more reasonable amplitudes than the SOSC pole figures. Finally,Fig. 16(d) shows the texture predicted by D-Rex with recrystallization. The orientations of the [100] and [010] axes now agree better with the laboratory data, but the amplitudes of the pole figures are too high by a factor of ≈2.

To summarize, in simple shear there is a range of strains for which the CPO is not well predicted by the finite strain hypothesis. This is an unavoidable limitation of the models presented in this paper, which characterize all deformation paths solely by the FSEs that they produce. However, we judge that this limitation is acceptable for two reasons. First, the finite strain hypothesis only fails for an intermediate range of strains, roughly γ = 50-300 per cent. Second, within this range the maximum departure of the [100] axis orientation from the long axis of the FSE is only about 25 • [START_REF] Zhang | Lattice preferred orientation of olivine aggregates deformed in simple shear[END_REF], which should be compared with the maximum possible error of 90 • .

A second limitation of our models is their inability to model textural evolution by reactivation of a pre-existing CPO, or (more generally) when the macrosopic velocity gradient tensor experienced by the aggregate changes rapidly with time. On the basis of laboratory experiments and observations from natural shear zones, [START_REF] Skemer | The influence of deformation history on the interpretation of seismic anisotropy[END_REF] concluded that pre-existing CPO introduces a lag between a change in the deformation kinematics and realignment of the olivine [100]-axis with the shear direction. [START_REF] Boneh | The effect of deformation history on the evolution of olivine CPO[END_REF] deformed dunite samples with strong pre-existing textures in three different configurations, with the pre-existing foliation perpendicular, oblique and parallel to the shortening axis of the deformation. They found that the new CPO depended strongly on the orientation of the preexisting texture, and evolved differently in each experiment. Clearly a finite-strain parametrization is unable to describe such behaviour. However, the cases studied by [START_REF] Boneh | The effect of deformation history on the evolution of olivine CPO[END_REF] correspond to an essentially instantaneous change of deformation type, and as such represent an extreme situation that is not likely to occur in smoothly varying flow fields like those of mantle convection. But even if not instantaneous, deformation changes in mantle flow fields may still be rapid in certain regions, such as where the flow 'turns the corner' in classical corner-flow models of ridges and subduction zones [START_REF] Mckenzie | Finite deformation during fluid flow[END_REF]. In such conditions the finite strain hypothesis is likely to break down locally [START_REF] Blackman | Seismic anisotropy in the upper mantle 2. Predictions for current plate boundary flow models[END_REF], and so finite strain-based predictions of CPO in such regions should be interpreted with caution.

A final limitation of our models is that they do not incorporate the influence of orthopyroxene (enstatite), a secondary phase whose presence tends to reduce the overall anisotropy of deforming olivine-dominated rocks [START_REF] Kaminski | D-Rex, a program for calculation of seismic anisotropy in the convective upper mantle[END_REF]. Inclusion of enstatite in our models will be the subject of future work.

In conclusion, we suggest that the methods developed in this study may prove useful in modelling texture evolution in other kinematically deficient minerals that have too few slip systems to accommodate an arbitrary imposed deformation. An example is the lower-mantle phase MgSiO 3 post-perovskite, which has only two dominant slip systems [100](010) and [001](010) in addition to twinning modes that may help accommodate strain parallel to [100] and [010] [START_REF] Tommasi | Deformation, crystal preferred orientations, and seismic anisotropy in the Earth's D layer[END_REF]. Future work on this problem will employ a modified form of the evolution eq. ( 3) that contains a source term to model the effect of mechanical twinning [START_REF] Barton | The use of discrete harmonics in direct multi-scale embedding of polycrystal plasticity[END_REF].

Figure 1 .

 1 Figure 1. Two conventions for defining the Euler angles. (a) Bunge's convention for the angles (ϕ 1 , , ϕ 2 ). (b) Slip-system specific Euler angles (φ, θ , ψ),where l is the slip direction and n is the normal vector to the slip plane. x, y and z are the coordinate axes in the frame of the sample.

Figure 2 .

 2 Figure 2. Pole figures of the SBFs F 1 (top), F 2 (middle) and F 3 (bottom) for uniaxial compression with enhanced strains R 12 = 0 and R 23 = -R 31 = 1.5.

  pole figures all have circular symmetry because the orientation of the compressional axis is at the centre of each figure.

  .1, 0.5, 1, 2}. The resulting SBFs are plotted in Fig. 4 along a line in the Euler space with ϕ 1 = 0, = 90 • and ϕ 2 ∈ [90 • , 270 • ]. The increasing sharpness of the SBF with increasing deformation is evident.

Figure 3 .

 3 Figure 3. Pole figures of the SBFs F 1 (top), F 2 (middle) and F 3 (bottom) for plane strain with enhanced strains R (s) 12 = R (s) 23 = and R (s) 31 = -2. This deformation consists of extension along the x-axis and an equal amount of compression along the z-axis.

FFigure 4 .

 4 Figure 4. Behaviour of the SBF F 1 as a function of increasing strain, for uniaxial compression with enhanced finite strain parameters R (1) 12 = 0 and R (1) 23 = {0.1, 0.5, 1, 2}. The SBFs are shown as functions of ϕ 2 along a line ϕ 1 = 0, = 90 • in the Euler space.

Figure 5 .

 5 Figure 5. Pole figures of the difference ODF F 1 -F 2 -F 3 for uniaxial compression with p 12 = p 23 = 0, R 12 = 0 and R 23 = 0.15. The pole figures are to be compared with (24).

Figure 6 .

 6 Figure 6. Instantaneous slip rate γ 1 (θ , ψ) for the slip system (010)[100], predicted by the SOSC model with p 12 = p 23 = 0 for an initially isotropic olivine aggregate deformed in uniaxial compression along the z-axis with a shortening rate ˙ 0 ≡ -˙ 3 > 0. (a) ˙ 0 t = 0; (c) ˙ 0 t = 0.4. Panels (b) and (d) show the function (28) for A 1 = 1.66 and 2.28, respectively. Colour scale is in units of ˙ 0 .

Figure 7 .

 7 Figure 7.Amplitudes A 1 , A 2 and A 3 of the single-slip rotation rates that best fit the predictions of the SOSC model with slip systems of equal strength (p 12 = p 23 = 0), as functions of the equivalent strain r 0 . Results are shown for deformation in plane strain (˙ 1 = -˙ 3 = ˙ 0 , ˙ 2 = 0) and uniaxial compression (˙ 3 = -˙ 0 , ˙ 1 = ˙ 2 = ˙ 0 /2).

Figure 8 .

 8 Figure 8. Amplitude A 1 = A 2 = A 3 of the single-slip rotation rate that best fits the predictions of the SOSC model with equal slip system strengths (p 12 = p 23 = 0), as functions of the logarithmic axial ratios r 12 and r 23 of the finite strain ellipsoid. (a) Best-fitting amplitude of the single-slip model. (b) Analytical fit to the function shown in part (a) (Appendix D).

  Figs 9(a)-(c) show the amplitude factors H s ( p 12 , p 23 ) that give the best fit to the SOSC predictions. Again, in order to use these results we must fit the three factors H s to analytical functions of p 12 and p 23 . The explicit forms of the functions we use are given in Appendix D, and contour plots of them are shown in Figs 9(d)-(f).

Figure 9 .

 9 Figure 9. Amplitude factors H s (p 12 , p 23 ) that measure the effect of unequal slip system strengths on the rate of crystallographic rotation for an undeformed (isotropic) aggregate with r 12 = r 23 = 0. Parts (a)-(c) show the amplitude factors predicted by the SOSC model for slip systems s = 1, 2 and 3, respectively. Parts (d)-(f) show the corresponding analytical fits to the data of parts (a)-(c) (Appendix D).

Figure 10 .

 10 Figure 10. (a) Function K defined by (31), obtained by numerical integration of the amplitude shown in Fig. 8(b). (b) Bicubic polynomial fit (D8) to the function of part (a).

Figure 11 .

 11 Figure 11. Expansion coefficients C 1 , C 2 and C 3 (solid lines) as a function of equivalent strain r 0 for slip model 1 and deformation in (a) uniaxial compression and (b) plane strain. The CRSS ratios are τ 1 /τ 2 = τ 2 /τ 3 = 1 (p 12 = p 23 = 0) in both cases. The dashed lines show the variance reduction VR (east scale) of the fit of the expansion (36) to the numerical solution of the evolution eq. (3). The short vertical black bars indicate the critical strain above which the ODF expansion (36) starts to take on negative values in portions of the Euler space.

Figure 12 .

 12 Figure 12. Comparison of pole figures of (a) the ODF determined by the method of characteristics and (b) its best-fitting SBF expansion, for deformation in uniaxial compression to an equivalent strain r 0 = 0.27. The CRSS ratios are p 12 = p 23 = 0 in both cases, and the spin amplitudes A s = 1.663H s ( p 12 , p 23 )J (r 12 , r 23 ) are those of spin model 1. The compressional axis is in the centre of the diagram.

Figure 13 .

 13 Figure 13. Same as Fig. 12 but for plane strain with r 0 = 0.20.

Figure 14 .

 14 Figure 14. Pole figures for olivine aggregates deformed in uniaxial compression (compression axis horizontal). (a) Experimental data from Nicolas et al. (1973) for 58 per cent shortening (non-recrystallized grains only). The pole figures were created by approximating each orientation as a Gaussian distribution with half-width = 10 • . (b) Prediction of the SOSC model for 50 per cent shortening, p 12 = -0.693, p 23 = -0.405 and = 7.5 • . (c) Prediction of the SBF-2 model for 58 per cent shortening, p 12 = -0.693 and p 23 = -0.405. (d) Prediction of the D-Rex model (without recrystallization) for 58 per cent shortening, p 12 = -0.693, p 23 = -0.405 and = 7.5 • . Note the different scales for the colours in the different sets of pole figures.

  C s ( p 12 , p 23 , r 12 , r 23 ) H s ( p 12 , p 23 )K(r 12 , r 23 )r i j for spin model 1 or R (s) i j = H s ( p 12 , p 23 )r i j for spin model 2, using the explicit expressions for H s and K given in Appendix D.

  i j are defined by 15. (6) Rotate the SBFs according to the orientation of the FSE. (7) Calculate the final solution f (g) = 1 + 3 s=1 C s Fs (go) where Fs = F s -1.

Figure 16 .

 16 Figure 16. Pole figures corresponding to simple shear deformation to a shear strain γ = 140 per cent. The red lines indicate the shear plane, the dashed black lines show the orientation of the axis of maximum instantaneous extension, and the solid black lines show the orientation of the long axis of the finite strain ellipsoid. The sense of shear is shown by the arrows between parts (a) and (b). (a) Laboratory measurements of Lee et al. (2002). Pole figures were created by approximating each orientation as a Gaussian distribution with half-width = 10 • . (b) Predictions of the SOSC model with p 12 = -0.693, p 23 = -0.405 and = 7.5 • . (c) Predictions of the SBF-2 model with p 12 = -0.693 and p 23 = -0.405. (d) Predictions of the D-Rex model including recrystallization, with p 12 = -0.693, p 23 = -0.405, λ * = 5, M * = 125, χ = 0.2 and = 7.5 • .

Table 1 .

 1 Slip systems.

	Index s	Slip plane	Slip direction	τ s /τ 2	Exponent
	1	(010)	[100]	0.25-4.0	3.5
	2	(001)	[100]	1.0	3.5
	3	(010)	[001]	0.25-4.0	3.5
	4	(101)	[-101]	100.0	3.5
	5	( 1 0 -1)	[101]	100.0	3.5

Table 2 .

 2 Models used.

	Model	Numerical method	Spin
	SOSC	Self-consistent	Increases with strain
		homogenization
	CHAR-1	Characteristics	Model 1 (same as SOSC)
	CHAR-2	Characteristics	Model 2 (no increase with strain)
	SBF-1	SBF expansion	Model 1
	SBF-2	SBF expansion	Model 2

Table E1 .

 E1 Coefficients Ĉ(1)i jkl for the SBF-1 model.

	APPENDIX F: EXPANSION COEFFICIENTS FOR THE SBF-2 MODEL
	kl	ij = 00	ij = 01	ij = 10	ij = 02	ij = 11	ij = 20
		1.819	0.01951	0.03516	-0.008173	0.003625	0.01318
		-0.2144	-0.07953	-0.007785	0.02354	0.001792	-0.003952
		-0.2292	0.07315	-0.01130	-0.02110	-0.004380	-0.005089
		2.575	0.2239	0.1249	-0.05780	0.001298	0.04917
		8.268	0.2469	0.4707	-0.08008	0.04004	0.1645
		2.627	-0.1266	0.1295	0.05170	0.02958	0.05083
		-5.290	-0.4313	-0.2900	0.1187	-0.004003	-0.1066
		-11.83	-1.598	-0.8329	0.4483	-0.03716	-0.2982
		-11.92	0.7567	-0.8875	-0.2372	-0.1294	-0.3172
		-5.357	0.1726	-0.2969	-0.07471	-0.05759	-0.1098
		4.418	0.4176	0.2214	-0.1188	-0.006062	0.07866
		6.652	1.638	0.3021	-0.4337	-0.01883	0.1134
		42.12	1.629	3.084	-0.4999	0.2704	1.068
		6.677	-1.361	0.3636	0.4446	0.1095	0.1356
		4.438	-0.2101	0.2267	0.08251	0.05126	0.08117
	Table E2. Coefficients	Ĉ(2) i jkl for the SBF-1 model.
	kl	ij = 00	ij = 01	ij = 10	ij = 02	ij = 11	ij = 20
		0.1800	-0.01872	0.02625	0.007899	0.005443	0.007778
		-0.3897	0.05687	-0.01316	-0.01755	-0.008754	-0.002546
		0.3041	-0.04911	0.01619	0.01265	0.005439	0.004547
		0.2929	-0.001004	0.05517	0.003216	0.01405	0.009504
		-1.437	-0.09096	0.1848	0.04176	0.05198	0.04708
		3.299	-0.2216	-0.3236	0.06998	0.01102	-0.09897
		1.4156	-0.01793	-0.1202	-0.01309	-0.04691	-0.03432
		-4.782	0.6631	-0.3276	-0.2210	-0.1526	-0.05571
		2.370	-0.5683	-0.2584	0.1740	-0.005765	-0.04158
		-6.522	0.5918	0.8870	-0.1984	-0.02436	0.2544
		-1.678	-0.05997	0.08054	0.03583	0.05266	0.02709
		10.07	0.1489	0.1043	0.00005165	-0.004820	-0.01680
		-19.89	-1.309	1.733	0.3979	0.5015	0.4964
		5.619	0.6033	-0.2083	-0.1955	-0.05099	-0.1129
		6.809	-0.5115	-0.8944	0.1735	0.01113	-0.2570
	Table E3. Coefficients	Ĉ(3) i jkl for the SBF-1 model.
	kl	ij = 00	ij = 01	ij = 10	ij = 02	ij = 11	ij = 20
		0.1800	0.04297	0.02723	-0.01354	-0.004171	0.008085
		0.2910	0.07013	0.01806	-0.02439	-0.005761	0.005476
		-0.3839	-0.07454	-0.01133	0.02751	0.009842	-0.002128
		3.334	-0.05532	-0.3217	0.003076	-0.02746	-0.09954
		-1.431	0.2938	0.1938	-0.1228	-0.06237	0.04838
		0.3350	0.05093	0.04114	-0.02648	-0.01814	0.005145
		-6.560	0.2331	0.88412	-0.07114	0.03229	0.2550
		2.316	0.1747	-0.3082	0.04193	0.09141	-0.04789
		-4.817	-1.101	-0.3289	0.4411	0.2189	-0.05466
		1.237	-0.04551	-0.09107	0.02660	0.02566	-0.02515
		6.810	-0.3289	-0.8870	0.1071	-0.02132	-0.2559
		5.561	-0.5921	-0.1370	0.08768	-0.08971	-0.1007
		-19.83	3.321	1.898	-1.365	-0.5194	0.5412
		10.25	-0.08463	0.05676	0.04518	-0.06226	-0.03257
		-1.490	0.08246	0.06357	-0.03757	-0.03056	0.02174

Table F1 .

 F1 Coefficients Ĉ(1)i jkl for the SBF-2 model.

	kl	ij = 00	ij = 01	ij = 10	ij = 02	ij = 11	ij = 20
		1.819	0.01944	0.03519	-0.008150	0.003628	0.01319
		-0.02042	-0.02953	0.002288	0.008983	0.0008590	-0.00002744
		-0.03022	0.03415	0.0007045	-0.009207	-0.0009204	-0.0004993
		0.4352	-0.007944	0.01594	0.003808	0.002214	0.007931
		1.798	0.03632	0.07473	-0.01536	0.005362	0.02710
		0.4589	0.01072	0.01498	-0.001385	0.002416	0.007483
		-0.4674	0.04018	-0.02626	-0.01130	-0.004565	-0.01082
		0.02793	-0.1278	0.01301	0.03685	0.004658	0.0009552
		-0.005665	0.1588	0.005230	-0.04495	-0.003242	-0.001527
		-0.4897	-0.05731	-0.02390	0.01723	-0.0002785	-0.009969
		0.2418	-0.01124	0.01346	0.002585	0.001585	0.005039
		-1.077	0.04456	-0.07192	-0.008495	-0.007175	-0.02236
		2.763	0.09402	0.1762	-0.03528	0.01335	0.06098
		-1.061	-0.1337	-0.06748	0.04239	-0.002232	-0.02092
		0.2484	0.02289	0.01248	-0.007143	0.0005071	0.004700
	Table F2. Coefficients	Ĉ(2) i jkl for the SBF-2 model.
	kl	ij = 00	ij = 01	ij = 10	ij = 02	ij = 11	ij = 20
		0.1804	-0.01867	0.02621	0.007884	0.005437	0.007766
		-0.2798	0.03601	0.001183	-0.01090	-0.004048	0.001009
		0.3038	-0.03864	-0.006664	0.01044	0.003672	-0.002039
		0.2133	0.0002898	-0.005922	-0.0004338	-0.001107	-0.004053
		-0.2839	-0.03413	0.005726	0.01616	0.008572	0.001882
		0.4155	-0.02899	-0.006770	0.008756	0.003378	-0.004159
		0.2384	-0.02207	-0.001968	0.00485	0.0008164	-0.00002040
		-1.732	0.1903	0.06083	-0.05739	-0.01699	0.02235
		0.4241	-0.1775	0.02248	0.04969	0.02492	0.008219
		-0.05844	0.04923	-0.005457	-0.01655	-0.006589	-0.001346
		-0.2093	0.007717	0.003461	-0.0009428	0.0008611	0.001284
		1.684	-0.04852	-0.05471	0.01218	-0.003299	-0.02070
		-1.706	-0.08726	0.08337	0.03327	0.02769	0.02778
		0.4752	0.1166	-0.03685	-0.03693	-0.02042	-0.01433
		0.05517	-0.02422	-0.006496	0.008482	0.002714	-0.001689
	Table F3. Coefficients	Ĉ(3) i jkl for the SBF-2 model.
	kl	ij = 00	ij = 01	ij = 10	ij = 02	ij = 11	ij = 20
		0.1802	0.04305	0.02719	-0.01357	-0.004180	0.008073
		0.2931	0.03646	-0.004282	-0.01346	-0.004464	-0.001152
		-0.2705	-0.03982	0.001700	0.01449	0.004989	0.001108
		0.4402	0.02964	-0.01120	-0.01322	-0.007001	-0.005680
		-0.2815	0.03701	0.0009341	-0.01183	-0.005890	0.0003152
		0.1947	0.001551	-0.007817	-0.002628	-0.002620	-0.004662
		-0.08161	-0.05455	0.0002216	0.01789	0.005398	0.0003776
		0.3844	0.2135	0.03621	-0.07683	-0.02383	0.01325
		-1.704	-0.1570	0.06717	0.05992	0.02528	0.02401
		0.2547	0.01906	-0.0003668	-0.006217	-0.001404	0.0005521
		0.06197	0.01689	-0.008338	-0.004968	-0.001841	-0.002234
		0.4937	-0.1499	-0.04303	0.04954	0.01274	-0.01679
		-1.704	0.1654	0.08456	-0.05811	-0.01755	0.02816
		1.674	0.009837	-0.05950	-0.006101	-0.006163	-0.02211
		-0.2125	-0.006105	0.003143	0.002333	0.0003854	0.001154

A C K N O W L E D G E M E N T S

We thank T. Becker, C. Conrad, A. Davaille, M. Faccenda, G. Ito, M. Long, L. Mameri, J. Park and A. Tommasi for helpful discussions. A. Tommasi kindly provided a Matlab script for plotting pole figures, as well as the crystal orientation data used to make Fig. 14(a). S.-I. Karato and H. Jung generously sent us the data used to make Fig. 16. The comments and suggestions of two anonymous referees helped us greatly to improve the original manuscript. This work was supported by grant PTECTO from the French Agence Nationale de la Recherche (ANR) and by the Programme National de Planétologie (PNP) of the Institut National des Sciences de l'Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS), co-funded by the Centre National d'Etudes Spatiales (CNES). NMR is grateful for the hospitality of the Isaac Newton Institute at the University of Cambridge where part of this work was done, and thanks the Simons Foundation for financial support during that time.

A P P E N D I X A : A N A LY T I C A L D E TA I L S O F T H E S I N G L E -A N D T R I P L E -S L I P M O D E L S A1 Single-slip model

The slip rate γ is determined by minimizing the misfit E defined by (6), either directly (unconstrained minimization) or subject to the strain rate compatibility constraint (8). Consider first the unconstrained case. Defining L = (e i j -E i j )(e i j -E i j ), e i j = γ S i j , (A1)

we can write down the Euler equation of the calculus of variations as

which implies

Solving (A3) for γ and noting that S ij S ij = 1/2, we obtain γ = 2S i j E i j .

(A4)

Turning now to the constrained case, the Euler equation analogous to (A2) is

where λ ij are Lagrange multipliers and λ ji = λ ij . The strain rate compatibility constraint is

Evaluating the derivatives in (A5) and solving for γ , we obtain

The integrals K ijkl can be evaluated explicitly by setting l i = a i1 and n i = a 2i , where a ij are the direction cosines (1). Eq. (A8) comprises five simultaneous equations for the five independent Lagrange multipliers λ ij , and has the solution λ ij = -4E ij . Substituting this result into (A7), we obtain γ = 10S i j E i j .

(A9)

A2 Triple-slip model

Consider an olivine aggregate each crystal of which has the three slip systems (010)[100] (s = 1), (001)[100] (s = 2), and (010)[001] (s = 3). Suppose further that all three slip systems have the same critical resolved shear stress (p 12 = p 23 = 0). As for the single-slip case, we begin with unconstrained minimization. The analogue of (A1) is

where γs is the slip rate on slip system s and

is the Schmid tensor for slip system s, defined by unit vectors l (s) in the slip direction and n (s) normal to the slip plane. The components of these vectors are

where a ij are the direction cosines (1). The Euler equations are

The solutions of these equations are γs = 2S

(s)

which is formally identical for each slip system to the unconstrained single-slip model prediction (A4). However, the predictions of the constrained single-slip and triple-slip models are very different. The analogue of (A5) is

the solutions of which are

Evaluating the integrals K (s)

i jkl explicitly with the help of (A11) and (A12), we find that (A18) has the solution

Comparing the preceding result with the canonical form γs = 2A s S (s)

i j E i j (see eq. ( 9)), we find A s = 5/3. This agrees to within 0.2 per cent with the value A s = 1.663 that we found by fitting the expression γs = 2A s S (s)

i j E i j to the SOSC model predictions for the case p 12 = p 23 = 0.

A P P E N D I X B : A N A LY T I C A L S O L U T I O N O F T H E S I N G L E -S L I P E V O L U T I O N E Q UAT I O N

The first step is to note that ( 14) implies ∂/∂t = A∂/∂T, whence (17) becomes

The characteristics satisfy

Integrating (B2) subject to the condition ψ(0) = ψ 0 , we obtain

where

Now the solution along the characteristics is

Now use the expression (B3a) for tan ψ to rewrite cos 2ψ and sin 2ψ in the expression for ∂ ψ/∂ψ. The integral in (B4) can then be evaluated to yield

Finally, use (B3a) to eliminate tanh -1 from (B5), and substitute the result into the expression (B4) for f. After some rearrangement, we obtain (18).

A P P E N D I X C : C RY S TA L L O G R A P H I C RO TAT I O N R AT E S

The explicit expressions for the rotation rates ġs ≡ ( φs , ψs , θs ) are

where

)

)

)

)

)

)

and (φ, ψ, θ) are the slip system-specific Eulerian angles (Fig. 1b) for slip system (010)[100] (s = 1).

A P P E N D I X D : A M P L I T U D E O F T H E C RY S TA L L O G R A P H I C RO TAT I O N R AT E

In this appendix we quantify the dependence of the slip system amplitudes A s (p 12 , p 23 , r 12 , r 23 ) on the CRSS ratios p 12 and p 23 and the strains r 12 and r 23 .

D1

Limit A s (0, 0, r 12 , r 23 )

Consider first how the activities of slip systems with equal strengths (p 12 = p 23 = 0) vary as a function of strain. In this case, the SOSC model shows that the amplitudes are equal for all slip systems, that is,

where the factor 1.663 is the amplitude when r 12 = r 23 = 0. Using a least-squares procedure, we find that the data of Fig. 8(a) are well fit if J (x) = exp 0.2891x 2 + 12.41x 3 -31.14x 4 + 27.06x 5 , (D2a) We verified that the above results are independent of deformation type by performing additional calculations for plane strain (˙ 1 = -˙ 3 , ˙ 2 = 0), and found that the fitted amplitudes A 1 -A 3 differ by < 0.1 per cent from those determined for uniaxial compression.

D3 General case

A general form for A s (p 12 , p 23 , r 12 , r 23 ) that is consistent with both limiting results (D2) and (D6) is