Dénes Takács
email: takacs@mm.bme.hu

Gábor Orosz

Gábor Stépán
email: stepan@mm.bme.hu

Delay effects in shimmy dynamics of wheels with stretched string-like tyres

Keywords: shimmy, delay, relaxation length, travelling wave solution, Hopf bifurcation

The dynamics of wheel shimmy is studied when the self-excited vibrations are related to the elasticity of the tyre. The tyre is described by a classical stretched string model, so the tyre-ground contact patch is approximated by a contact line. The lateral deformation of this line is given via a non-holonomic constraint, namely, the contact points stick to the ground, i.e., they have zero velocities. The mathematical form of this constraint is a partial differential equation (PDE) with boundary conditions provided by the relaxation of deformation outside the contact region. This PDE is coupled to an integro-differential equation (IDE), which governs the lateral motion of the wheel. Although the conventional stationary creep force idea is not used here, the coupled PDE-IDE system can still be handled analytically. It can be rewritten as a delay differential equation (DDE) by assuming travelling wave solutions for the deformation of the contact line. This DDE expresses the intrinsic memory effect of the elastic tyre. The linear stability charts and the corresponding numerical simulations of the nonlinear system reveal periodic and quasi-periodic selfexcited oscillations that are also confirmed by simple laboratory experiments. The observed quasi-periodic vibrations cannot be explained in single degree-of-freedom wheel models subject to a creep force.

Introduction

The lateral vibration of towed wheels, called shimmy, may appear on airplane landing gears, motorcycle wheels, caravans, rear wheels of semi-trailers and articulated buses, and it usually presents a safety hazard. The terminology 'shimmy' appeared in the 1920's when it was the name of a popular dance. The shimmy of towed wheels may be caused either by the elasticities of the towing bar suspension and the attached vehicle structure, or by the elasticity of the tyre on the wheel, or by a combination of the two cases. There are two important reasons why the shimmy phenomenon has not been fully explored yet. One part of the problem is that the vehicle itself is a complex dynamical system serving several low-frequency vibration modes which may all be important components of the dynamical behaviour at different running speeds and conditions. This part of the problem is mainly resolved by the development of multibody dynamics and commercial computer codes [START_REF] Schiehlen | Computational dynamics: theory and applications of multibody systems[END_REF]. The other part of the problem is originated in the wheel-ground contact. Advanced finite element calculations need long computation times even in stationary cases; see [START_REF] Kalker | Wheel rail rolling-contact theory[END_REF] for railway wheels and [START_REF] Böhm | Model for the radial tire for high frequent rolling-contact[END_REF][START_REF] Chang | Literature survey of transient dynamic response tyre models[END_REF] for tyres. Even nowadays, the dynamic contact problems usually require special codes, large computational power, and still, there are no analytical results available to check these calculations.

One of the first scientific reports on shimmy was presented by von Schlippe and Dietrich [START_REF] Schlippe | Dynamics of flexible multibody systems having rolling contact: application of the wheel element to the dynamics of road vehicles[END_REF], where they analysed a simple low degree-of-freedom (DoF) model. In this so-called stretched string model the tyre-ground contact was considered as a contact line that becomes deformed due to the lateral displacement of the wheel (see the curve between the points L and R in Fig. 1), while the longitudinal deformation was neglected. Furthermore, it was considered that each contact point P sticks to the ground. Note that the tyre becomes deformed not only along the contact line but also in front of the leading point L and behind the rear point R as represented by the deformed central line of the tyre in Fig. 1. It was assumed that outside the contact patch the deformation decays exponentially, which was also confirmed by measurements. Because the resulting equations were too complicated to be analysed with the available mathematical tools at that time, they introduced a severe simplification, namely, the contact line was straight between the points L and R. This way, the resultant lateral force and torque induced by the elastic tyre deformation were calculated, leading to a delay differential equation (DDE) with a discrete delay. Using this equation the linear stability of the stationary rolling motion was analysed (with some further simplifications since the mathematical theory of DDEs was not available at that time). The discrete delay was equal to the time period of a tyre point while in contact with the ground between L and R. Since then, several versions of the stretched string model have been developed and analysed. For example, in [START_REF] Segel | Force and moment response of pneumatic tires to lateral motion inputs[END_REF] A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT frequency response functions were calculated for the stretched string model without any restriction to the shape of the contact line. This so-called exact stretched string model has become a basic reference for later studies. This approach results in a DDE with a distributed delay as explained in details in Section 3.

A different approach was taken by Pacejka [START_REF] Pacejka | The wheel shimmy phenomenon[END_REF] who introduced different straight and curved contact line approximations by using stationary shape functions for the lateral tyre deformation calculated at constant drift angles. The resultant lateral creep force and torque are calculated from these 'quasi-stationary' deformations by the semi-empirical 'Magic Formula' [START_REF] Pacejka | Tyre and Vehicle Dynamics[END_REF]. This way, the delay effects are completely eliminated. In the resulting simplified models only the caster angle and the lateral deformation of the leading point L were used as state variables, leading to an ordinary differential equation (ODE) that made it very popular and easy-to-analyse. In the middle range of the towing speeds, the quasi-stationary deformation idea gives reasonable agreement with experiments. Several research reports prove the success of this approach in engineering (see [START_REF] Troger | A nonlinear-analysis of the generic types of loss of stability of the steady-state motion of a tractor-semitrailer[END_REF] on tractor-semi-trailers systems, [START_REF] Sharp | Advances in the modelling of motorcycle dynamics[END_REF] on motorcycles or [START_REF] Fratila | Simulation of coupled car and caravan handling behaviour[END_REF] on caravans). Nonlinearities were also introduced and their importance were emphasized by the existence of unstable periodic motions in [START_REF] Pacejka | The wheel shimmy phenomenon[END_REF]. The model was also generalised by considering the effect of the width of the contact area, that of sliding at the rear part of the contact patch, and that of the gyroscopic effects appearing when the wheel is allowed to be tilted from its vertical plane. Pacejka's wisdom about tyres has accumulated in his book [START_REF] Pacejka | Tyre and Vehicle Dynamics[END_REF] that also includes a separate chapter for shimmy with an extensive reference list.

Considering the exact stretched string model with simplified boundary conditions, Stépán [START_REF] Stépán | Delay, nonlinear oscillations and shimmying wheels[END_REF] has introduced nonlinearities into the system. He also investigated the linear stability with mathematical rigor leading to the possibilities of quasi-periodic oscillations, which recently has been confirmed experimentally in [START_REF] Takács | Dynamics of Rolling of Elastic Wheels[END_REF]Takács and Stépán, to appear 2007).

The elasticity of the tyre was also considered in a point contact model by Moreland [START_REF] Moreland | The story of shimmy[END_REF]. The contact line was shrunk into a point where the force, induced by the elastic tyre, acts. A relaxation time was also introduced for the force to model its 'delayed action' and a torque coefficient was defined to relate the force and the torque. It was proven by Collins [START_REF] Collins | Theories on mechanics of tyres and their application to shimmy analysis[END_REF] that the point contact model is equivalent to the stretched string model when the latter is restricted to the case of straight contact line. However, in the point contact model the relaxation time and a torque coefficient has to be estimated or measured, while the stretched string model provides the corresponding constants via its geometry.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

If there is elasticity in the suspension system, even a point contact model with rigid tyre can exhibit shimmy and complicated nonlinear (sometimes even chaotic) behaviour [START_REF] Goodwine | Controlling unstable rolling phenomena[END_REF][START_REF] Saux | Dynamics of a rolling disk in the presence of dry friction[END_REF][START_REF] Schlippe | Dynamics of flexible multibody systems having rolling contact: application of the wheel element to the dynamics of road vehicles[END_REF][START_REF] Stépán | Chaotic motion of wheels[END_REF][START_REF] Stépán | Appell-Gibbs equation for classical wheel shimmy -an energy view[END_REF]Takács et al., 2007). However, feed-back linearisation based controllers can be constructed for these nonlinear systems to suppress the vibrations [START_REF] Goodwine | Feedback stabilization of a class of unstable nonholonomic systems[END_REF].

In this study, we consider the case when the wheel is pulled by a caster fixed to a cart of constant velocity. The lateral deformations of the tyre are modelled by the exact stretched string model. The tyre-ground contact is described by a contact line. We assume that each contact point sticks to the ground that results in a non-holonomic constraint expressed by a first-order partial differential equation (PDE). We also assume that the deformation decays exponentially outside the contact region with a characteristic relaxation length. With the appropriate choice of the boundary conditions, the relaxation length of the tyre is taken into account among other conventional tyre parameters like the specific stiffness and damping. The Newtonian equation of motion becomes a second order integro differential equation (IDE). Assuming travelling wave solutions of the deformation allows us to transform the PDE-IDE system into a delay differential equation (DDE) with distributed delay. Here, the delay is the time needed for the leading point L of the contact line to travel backward (relative to the caster) to the actual contact point P (see Fig. 1). The linear stability investigation of the DDE shows that the stationary rolling motion may lose its stability via co-dimension one or co-dimension two Hopf bifurcations as the parameters (like the towing speed and the caster length) are varied. Consequently, self-excited periodic and quasi-periodic oscillations can appear. The stability chart in the plane of the above parameters is determined analytically and checked by numerical simulations and laboratory experiments.

Mechanical model

Consider the simple model of a towed wheel with elastic tyre in Fig. 1. The wheel is pulled by a caster, and the suspension point A of the caster is towed with a constant velocity v. The length of the caster is l and the tyre-ground contact length is 2a. We consider the suspension system to be rigid and use the exact stretched string model of the tyre.

One of the chosen state variables is clearly the caster angle ψ of rotation about the vertical axis. In accordance with the stretched string model, the tyre-ground contact patch is approximated by the contact line, i.e., by the deformed central line of the tyre between the leading point L and the rear point R (see Fig. 1). Note that the deformation of the tyre outside the contact line is

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT (a) X Y Z x y z ψ(t) A R P L v (b) x z v A R P L X Y x y ψ(t) q(x, t) v l a a A R P L
Fig. 1. Model of a towed wheel with elastic tyre. Panel (a) shows the 3-dimensional view of the wheel while panel (b) depicts the side and top view of the wheel. The deformed central line of the tyre is shown in both panels. This forms the contact line between the points L and R where the tyre is connected to the ground. The (x, y, z) coordinate system is fixed to the caster while the (X, Y, Z) coordinate system is fixed to the ground.

also described by the deformation of its central line. The lateral deformation q(x, .) of the central line is described in the coordinate system (x, y, z) fixed to the caster as shown in Figs. 1 and2. This is a state variable distributed along the contact line x ∈ [-a, +a], while it is defined by the exponentially decaying functions

q(x, t) = ⎧ ⎨ ⎩ q(a, t) e -(x-a)/σ , if x ∈ [a, ∞) , q(-a, t) e (x+a)/σ , if x ∈ (-∞, -a] , (1)
before the leading point L and behind the rear point R. The tyre parameter σ is called the relaxation length (see Fig. 2(b)). This is considered to be small relative to the wheel diameter, so the theoretical extension x ∈ (-∞, +∞) is an acceptable standard approximation [START_REF] Pacejka | Tyre and Vehicle Dynamics[END_REF][START_REF] Schlippe | Dynamics of flexible multibody systems having rolling contact: application of the wheel element to the dynamics of road vehicles[END_REF][START_REF] Segel | Force and moment response of pneumatic tires to lateral motion inputs[END_REF]. For general theory of deformation of strings see [START_REF] Kármán | Mathematical Methods in Engineering[END_REF], and for derivation of the equations of motion for moving continua see [START_REF] Wickert | Classical vibration analysis of axially moving continua[END_REF].

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT (a) X Y Z x y z L P R (b) x y +a -a +(a + σ) -(a + σ) R L P q(x, t)
Rolling without sliding means a non-holonomic constraint with respect to the state variables ψ and q(x, .) which is formulated as follows. In the ground-fixed coordinate system (X, Y, Z), the position vector of a contact point P is given as

⎡ ⎢ ⎣ X(x, t) Y (x, t) ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ vt -(l -x) cos ψ(t) -q(x, t) sin ψ(t) -(l -x) sin ψ(t) + q(x, t) cos ψ(t) ⎤ ⎥ ⎦ (2) for x ∈ [-a, a].
(The trivial condition Z(x, t) ≡ 0 is not spelled out here since the vertical dynamics are neglected.) Since point P sticks to the ground, its velocity col[ dX dt dY dt ] is zero, that yields

⎧ ⎪ ⎨ ⎪ ⎩ d dt q(x, t) = v sin ψ(t) + (l -x) ψ(t) , ẋ = -v cos ψ(t) + q(x, t) ψ(t) , (3) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT with d dt q(x, t) = q(x, t) + q (x, t) ẋ , (4) for x ∈ [-a, a].
Here d dt refers to total differentiation with respect to time t while dot and prime refer to partial differentiation with respect to time t and space x. By eliminating ẋ, we obtain the constraining equation for the state variables in the form of a first-order scalar partial differential equation (PDE)

q(x, t) = v sin ψ(t)+(l-x) ψ(t)+q (x, t) v cos ψ(t)-q(x, t) ψ(t) , x ∈ [-a, a] .
(5)

We assume that the deformation q(x, t) is continuously differentiable at the leading point L which is often referred to as 'no kink at L' condition [START_REF] Pacejka | The wheel shimmy phenomenon[END_REF]. This provides the mixed boundary condition

q (a, t) = - q(a, t) σ (6)
for the PDE (5) in accordance with the exponential decay of deformation in (1) in front of the leading point L.

The equation of motion of the wheel can be given as an integro-differential equation (IDE)

J A ψ(t) = -k ∞ -∞ (l -x)q(x, t)dx -b ∞ -∞ (l -x) d dt q(x, t)dx , ( 7
)
where J A is the mass moment of inertia of the wheel and the caster together with respect to the z axis at the king pin A, while k and b stand for the lateral stiffness and damping of the tyre per unit length, respectively. To simplify the model, the quantities k and b are considered to be constant even outside the contact region where the deformation of the wheel central line is projected to the ground as in Fig. 2 (even though these are constant rather along the circumference). Thus, the towed wheel model is given by the coupled PDE-IDE system (5,7) when considering the approximation (1) and the boundary condition (6).

Travelling wave solutions and time delay

Let τ (x) denote the time needed for a tyre particle located at the leading point L to travel backward (relative to the caster) to the actual point P characterized by the coordinate x. With this delay, we can describe a travelling wave solution

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT as ⎡ ⎢ ⎣ X(x, t) Y (x, t) ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ X(a, t -τ (x)) Y (a, t -τ (x)) ⎤ ⎥ ⎦ , ( 8
)
which simply expresses the fact that the tyre particles once fixed to the ground remain fixed there, while the caster travels ahead above them. With the help of (2), the travelling wave ( 8) can be transformed into the form

⎧ ⎨ ⎩ l -x = vτ cos ψ(t) + (l -a) cos ψ(t) -ψ(t -τ ) -q(a, t -τ ) sin ψ(t) -ψ(t -τ ) , q(x, t) = vτ sin ψ(t) + (l -a) sin ψ(t) -ψ(t -τ ) + q(a, t -τ ) cos ψ(t) -ψ(t -τ ) , ( 9
) which is still implicit due to the fact that τ (x) cannot be expressed in closed form. However, differentiating the first equation in ( 9) with respect to τ , we obtain

dx dτ = -v cos ψ(t) + (l -a) ψ(t -τ ) sin ψ(t) -ψ(t -τ ) -q(a, t -τ ) sin ψ(t) -ψ(t -τ ) +q(a, t -τ ) ψ(t -τ ) cos ψ(t) -ψ(t -τ ) . ( 10
)
Now, one may substitute (1,3,9) into ( 7) and use the change of variables x and τ (integration by substitution) based on (10). The resulting retarded functional differential equation (RFDE) contains the time dependent caster angle ψ(t) and its delayed values ψ(tτ ), and the time dependent leading point lateral deformation q(a, t) and its delayed values q(a, t-τ ). The constraining equation ( 5) with its boundary condition (6) provides an ordinary differential equation (ODE) for the leading point lateral deformation:

q(a, t) = v sin ψ(t) + (l -a) ψ(t) - v σ q(a, t) cos ψ(t) + 1 σ q 2 (a, t) ψ(t) . (11)
The above described nonlinear RFDE-ODE system was given explicitly for the special case of zero relaxation length (σ = 0) and zero damping (b = 0) in [START_REF] Stépán | Delay, nonlinear oscillations and shimmying wheels[END_REF]. For non-zero σ and b the explicit linearised equations are given in the next section.

Pacejka's creep force/moment model also uses the leading point lateral deformation as state variable, but instead of the travelling wave solution along the contact line, it uses the stationary lateral deformation obtained at constant drift angle ψ [START_REF] Pacejka | The wheel shimmy phenomenon[END_REF]. This way, one obtains a three dimensional nonlinear ODE instead of the infinite dimensional nonlinear RFDE-ODE, but loses the dynamics within the contact region, which may be important in certain parameter domains as identified later.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 4 Small oscillations around stationary rolling

The stationary rolling motion of the wheel is described by the trivial solution

ψ(t) ≡ 0 , q(x, t) ≡ 0 , x ∈ (-∞, +∞) . ( 12
)
Small shimmy oscillations around the stationary rolling can be described by the linearisation of the governing equations (7,9,10,11):

J A ψ(t) = -k a -a (l -x) q(x, t) + b k d dt q(x, t) dx -kσ(l -a -σ) q(a, t) + b k q(a, t) + v σ q(a, t) -kσ(l + a + σ) q(-a, t) + b k q(-a, t) - v σ q(-a, t) , ( 13
) ⎧ ⎨ ⎩ l -x = vτ + l -a ⇒ τ (x) = a-x v , q(x, t) = (vτ + l -a)ψ(t) -(l -a)ψ(t -τ ) + q(a, t -τ ) , ( 14
) dx dτ = -v , (15) q(a, t) = vψ(t) + (l -a) ψ(t) - v σ q(a, t) . ( 16
)
In ( 13) the integrals over the intervals (-∞, -a], and [a, ∞) are calculated in closed form by using the exponential decay functions in (1), while the remaining d dt q(x, t) can be obtained from the linearisation of (3). The second equation in ( 14) means that the lateral deformation q(x, t) can be expressed by the present and delayed values of the caster angle ψ and the leading point lateral deformation q(a, .).

In particular, considering the first equation in ( 14) at the rear point R, we obtain

x = -a ⇒ τ (-a) = 2a v . ( 17
)
and so the second equation in ( 14) at the rear point R gives

q(-a, t) = (l + a)ψ(t) -(l -a)ψ(t -2a v ) + q(a, t -2a v ) . ( 18
)
After substituting (14,15,17,18) into the IDE (13) and dividing the equation

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
with the mass moment of inertia J A , we obtain the form

ψ(t) + 2ζω n ψ(t) + ω 2 n ψ(t) = kv J A 2a v 0 (l -a + vτ )p(t -τ )dτ + kσ J A (l -a -σ) p(t) + 2ζ ω n ṗ(t) + v σ p(t) + kσ J A (l + a + σ) p(t -2a v ) + 2ζ ω n ṗ(t -2a v ) - v σ p(t -2a v ) + bv J A 2l(a + σ)ψ(t) , ( 19
)
where we temporarily introduced the new notation

p(t) = (l -a)ψ(t) -q(a, t) (20)
for the absolute position Y of the leading point L to shorten the expression. The constants

ω n = 2k J A a l 2 + a 2 /3 + σ l 2 + a 2 + aσ and ζ = ω n 2 b k (21)
are the undamped natural angular frequency and the damping ratio of the steady wheel (v = 0), respectively.

Rescaling

Let us rescale the time as

T := v 2a t , (22)
define the new integration variable by

ϑ := - v 2a τ , ( 23
)
and the dimensionless leading point lateral deformation by

Q(T ) := 1 a q(a, T ) . ( 24
)
Now, using dot as partial differentiation with respect to the rescaled time T , we can define the dimensionless angular velocity as

Ω(T ) := ψ(T ) . ( 25
) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
Furthermore, the dimensionless towing speed V , the dimensionless caster length L and the dimensionless relaxation length Σ are given by

V := 1 ω n v 2a , L := l a , Σ := σ a , (26) respectively.
Using definitions (22-26), the ODE ( 16) and the IDE (19,20) provide a 3dimensional system of first-order DDEs:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ψ(T ) Ω(T ) Q(T ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 0 -1 V 2 + c 1 c 2 -2ζ V -c 2 Σ 2 (L -1 -Σ) 2 L -1 -2 Σ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ψ(T ) Ω(T ) Q(T ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + c 2 0 -1 (L -1 -2ϑ) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 L -1 0 -1 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ψ(T + ϑ) Ω(T + ϑ) Q(T + ϑ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ dϑ + c 2 2 (L + 1 + Σ) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 (Σ-4ζV )(L-1)-4ΣζV 0 8ζV -Σ 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ψ(T -1) Ω(T -1) Q(T -1) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (27)
where

c 1 = Σ 2 (L -1 -Σ)(L -1) + 2ζV (L 2 + (1 + Σ) 2 ) , c 2 = 1 V 2 1 L 2 + 1/3 + Σ(L 2 + 1 + Σ) .
(28)

Stability investigation

The stationary rolling motion ( 12) is now represented by the trivial solution

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ψ(T ) Ω(T ) Q(T ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ≡ 0 (29) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
of the linearised equation of motion ( 27). The Laplace transformation of ( 27) or the substitution of the trial solution

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ψ(T ) Ω(T ) Q(T ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = Ke λT , K ∈ C 3 , λ ∈ C (30)
leads to the characteristic function

D(λ; μ) = ΣV 2 λ 3 + 2V (V + Σζ)λ 2 + (Σ + 4ζV )λ + 2 - L -1 -Σ L 2 + 1/3 + Σ(L 2 + 1 + Σ) × × 2 λ 2 (L -1)λ + 2 -(L + 1)λ + 2 e -λ + 4ζV L(1 + Σ)(2 + Σλ) L -1 -Σ + (L -1 -Σ)(2ΣζV λ + Σ + 4ζV ) + (L + 1 + Σ)(2ΣζV λ + Σ -4ζV )e -λ , (31)
where μ ∈ R 4 represents the dimensionless parameters ordered in a vector

μ = col [ V L Σ ζ ] . (32)
Generally, equation (31) has infinitely many complex zeros for the characteristic exponents (characteristic roots) λ, but only a finite number of these may be situated in the right-half complex plane. The stationary rolling (12,29) is asymptotically stable if and only if all the infinitely many characteristic exponents are situated in the left-half complex plane [START_REF] Stépán | Retarded Dynamical Systems: Stability and Characteristic Functions[END_REF].

At the limit of stability, bifurcation can take place in the corresponding nonlinear system when characteristic roots are located at the imaginary axis for some critical values μ cr of the parameter vector. It is easy to see, that the critical characteristic root cannot be the zero since it does not satisfy the characteristic equation for any parameter values, i.e., D(0, μ) = 0. This means that only Hopf bifurcation can occur at μ cr when a pair of pure imaginary complex conjugate characteristic exponents In the 4-dimensional parameter space μ ∈ R 4 , these formulae describe stability boundaries (3-dimensional hypersurfaces) parameterized by the dimensionless angular frequency ω ∈ R + . Fixing 2 of the 4 parameters, we obtain stability boundary curves in a parameter plane. In particular, we fix the dimensionless relaxation length Σ and the damping ratio ζ for different values, and construct stability boundary curves in the plane of the dimensionless towing speed V and the dimensionless caster length L. To decide whether a certain region bounded by the intricate structure of stability curves is stable or not, we use the stability criteria derived in [START_REF] Stépán | Retarded Dynamical Systems: Stability and Characteristic Functions[END_REF]. The stable parameter domains are shaded in the stability charts plotted in the (V, L)-plane in Fig. 3 for different values of Σ for the undamped system (ζ = 0). Notice that there exists a boundary at L = 1 + Σ for any Σ, and for large V the system is stable above this boundary and unstable below. Fig. 4 shows the stability charts in the (V, L)-plane for a mesured value of Σ (see Fig. 7 later) when ζ is increased. It can be observed that the unstable 'lenses' gradually shrink and disappear as the damping is increased and only the monotone increasing stability curve persist for large damping (ζ > 0.035). This curve saturates at L = 1 + Σ for large V .

λ 1,2 (μ cr ) = ±iω , ω ∈ R + (
Where the stability boundary curves intersect each other, two pairs of pure imaginary characteristic exponents ±iω 1 and ±iω 2 co-exist with two dimensionless angular frequencies ω 1 and ω 2 . This is a co-dimension two (or double) Hopf bifurcation which also arises in a similar delayed robot dynamics problem [START_REF] Stépán | Quasiperiodic oscillations in robot dynamics[END_REF]. Due to this bifurcation quasi-periodic self-excited oscillations appear around the stationary rolling motion with dimensional frequencies

f 1 = ω 1 v/(4aπ) = ω 1 V ω n /(2π) and f 2 = ω 2 v/(4aπ) = ω 2 V ω n /(2π).
The corresponding 'quasi-periodic travelling waves' in the contact region possess the dimensional wave lengths v/f 1 = 4aπ/ω 1 and v/f 2 = 4aπ/ω 2 .

Instead of carrying out the analytical study of co-dimension one and two Hopf bifurcations that would require the reduction of the dynamics from the infinite-dimensional state space to 2-and 4-dimensional centre manifolds, we identified typical periodic and quasi-periodic vibrations in the system by numerical simulations. These provide enough information at this stage of the research for the validation of our model by experiments.

V L V L V L V L (a) (b) (c) (d) ζ = 0 ζ = 0 Σ = 0 Σ = 0.5 ζ = 0 ζ = 0 Σ = 1.0 Σ = 1.8 L = 1 + Σ L = 1 + Σ L = 1 + Σ L = 1 + Σ

Simulation of the nonlinear equations

In order to demonstrate the stability properties determined above we study the original PDE-IDE system (5,7) with conditions (1,6) by numerical simulation.

On one hand we wish to verify that the obtained linear stability diagrams are correct. On the other hand we would like to obtain information about the appearing nonlinear oscillations when the stationary rolling motion is linearly unstable.

We fix the dimensional parameters a, σ, k and b as in Table 1 and vary v, l and J A such that ω n is kept constant (see ( 21)). Consequently, the dimensionless parameters Σ and ζ keep their values shown in Table 1 while V and L are varied (see ( 26)). This means that we consider the stability diagram in Fig. 4(c) where the points A-D are marked by crosses. We run the simulations using the parameters at these points. In each of the pairs A-C and B-D the points are separated by a Hopf curve such that one point lies in the stable (shaded) regime while the other in the unstable (white) regime. Consequently, qualitatively different behaviour is expected for the two points of each pair.

V L V L V L V L (a) (b) (c) (d) ζ = 0.005 ζ = 0.015 Σ = 1.8 Σ = 1.8 ζ = 0.02 ζ = 0.03 Σ = 1.8 Σ = 1.8 A C B D L = 1 + Σ L = 1 + Σ L = 1 + Σ L = 1 + Σ
To integrate the PDE-IDE system (5,7) we use the Lax-Wendroff method [START_REF] Lax | Systems of conservation laws[END_REF] which provides second order accuracy in time. Note that for the IDE component this method is effectively the same as the 2nd order Runge-Kutta method. The number of spatial mesh points was 400. Furthermore, we consider the initial condition

ψ(0) = 0 , ψ(0) = w , q(0, x) ≡ 0 , q(0, x) = (l -x)w , x ∈ [-a, a] , (35)
with constant w = 1[rad/s]. This corresponds to applying a lateral impact to the stationary rolling wheel at t = 0. In each panel the time profile for the caster angle ψ(t) is shown and the spatial distribution of the lateral deformation q * (x) = q(x, t * ) is depicted at time t * . In panels (a) and (c) the stationary rolling motion is linearly stable, while in panels (b) and (d) this motion is unstable and stable oscillating (shimmy) motions can be observed.

x[m] t[s] t * x[m] t[s] t * x[m] t[s] t * x[m] t[s] t * q * [m] ψ[rad] q * [m] ψ[rad] q * [m] ψ[rad] q * [m] ψ[rad] -π 2 π 2 -π 2 π 2 -π 2 π 2 -π 2 π 2 R L R L R L R L case A (a1) (a2) case B (b1) (b2) case C (c1) (c2) case D (d1) (d2)
The obtained results are shown in Fig. 5. In each panel the time history for the caster angle ψ(t) is shown together with the spatial distribution of the lateral deformation q * (x) = q(x, t * ) at the chosen time t * . Fig. 5(a) and (b) show the cases A and B where the stationary rolling motion is linearly stable, that is, small perturbations around this motion decay in time. Observe that the small amplitude oscillations/waves are very close to harmonic. Fig. 5(c) and (d) show the cases C and D where the stationary rolling motion is linearly unstable, that is, small perturbations grow in time and the system approaches largeamplitude oscillations as time progresses. In case C the approached motion is periodic corresponding to a limit cycle in phase space while in case D quasiperiodic oscillations are observed which correspond to a torus in phase space.

The large amplitude oscillations/waves are not harmonic anymore, but these are not realistic physically due to the large caster angles reaching even π/2. This will be further discussed in Section 8 on comparison to experimental observations. Table 1 The experimentally identified dimensional parameters and their dimensionless counterparts; see formulae ( 21) and ( 26). 2005; Takács and Stépán, to appear 2007). The wheel was placed on a conveyor belt of variable speed as shown in Fig. 6. In order to avoid oscillations due to the elasticity of the conveyor belt, it was stiffened laterally by a steel frame which also kept the possible lateral buckling of the conveyor belt under control. The caster length and the mass moment of inertia of the structure with respect to the vertical axis at the king pin A were also adjustable. In this way, we were able to tune the natural angular frequency ω n and the damping ratio ζ to desired values (see (21,26)). Thus, all the necessary parameters were controlled within certain limits to identify an experimental stability chart in the plane of dimensionless towing speed V and dimensionless caster length L (for fixed dimensionless relaxation length Σ and damping ratio ζ). Since all dimensionless parameters depend on the caster length l and the contact length 2a, the proper variation of the system parameters requires special attention.

In order to identify the numerical values of parameters, first, we fixed a certain air pressure in the pneumatic tyre, placed the wheel in a rigid frame as in Fig. 7(a) and pulled its centre point in lateral direction. enlarged (and distorted) picture of the deformed tyre in and around the contact region when a transparent plastic plate was placed at one side of the frame and the central line of the tyre was marked. This picture perfectly follows the approximation (1) of the stretched string model used in the literature [START_REF] Pacejka | Tyre and Vehicle Dynamics[END_REF][START_REF] Schlippe | Dynamics of flexible multibody systems having rolling contact: application of the wheel element to the dynamics of road vehicles[END_REF][START_REF] Segel | Force and moment response of pneumatic tires to lateral motion inputs[END_REF]. We identified the contact length 2a and the relaxation length σ, and the obtained results are shown in Table 1.

Then the standing wheel was placed to the conveyor belt as in Fig. 6 and its centre was slightly hit in lateral direction. The time history of the acceleration of a chosen caster point was recorded. From the frequency and logarithmic decrement of the vibration signal the natural angular frequency ω n and the damping ratio ζ can be determined. Using (21) the lateral stiffness k and lateral damping b per unit length can be calculated. The obtained results are also given in Table 1.

During the experiments, for a chosen caster length l, we increased the towing speed v step-by-step and identified the loss of stability of stationary rolling by detecting the appearance of self-excited vibrations, i.e., the shimmy. Then we repeated the same experiment for several different values of the caster length l. We managed to keep the constant value for the natural frequency f n when the caster length l was varied: we varied the mass at the end of the caster and the mass moment of inertia J A changed accordingly. Consequently, the damping ratio ζ and the dimensionless relaxation length Σ were also kept constant (see (21,26)).

The experimental stability chart was transformed to the (V, L)-plane of dimensionless parameters. the corresponding theoretical stability chart calculated from the time-delayed model for the experimentally identified parameters in Table 1 (see also the bottom part of Fig. 4(c)). This shows relatively good agreement between theory and experiment. Recall that the theoretical stability boundary saturates at L = 1 + Σ for large V . Fig. 8 shows that the experimental stability boundary saturates to a slightly higher value of L. This suggests that one may obtain a better fit for large V by considering slightly higher value of Σ than the measured one.

In Fig. 8 we marked by a cross the point D on the experimental stability boundary which is located close to the double Hopf point (the intersection of the theoretical stability curves). This is the same point as point D in Fig. 4(c) where quasi-periodic oscillations were found by numerical simulation as shown in Fig. 5(d). Our measurements confirm these observations as quasi-periodic vibrations appear in the experiments, too. However, we found that the amplitude of oscillations in the experiments is lower than suggested by simulations. These deviations are due to the shortcoming of our model that the tyre sticks to the ground even for very large lateral deformations, while in reality sliding usually occurs at the rear part of the contact region. Involving this dissipative effect one might be able to obtain better agreement between simulations and measurements. More details about this analysis can be found in (Takács and Stépán, to appear 2007).

Conclusion and discussions

A low degree-of-freedom model of the shimmying wheel with elastic tyre was investigated. The no-slip kinematic constraint along tyre-ground contact region was described by a nonlinear partial differential equation (PDE) which was coupled to an integro-differential equation (IDE) of wheel motion. Considering the relaxation of tyre deformation around the contact region provided a boundary condition for the coupled nonlinear PDE-IDE system.

The equations were linearized about the stationary rolling motion. Using travelling wave solutions the linearized PDE-IDE was transformed into a 3dimensional linear system of delay differential equations (DDEs). In this way, stability charts were constructed analytically in the plane of the towing speed and caster length for certain damping ratio and relaxation length parameters. Crossing the stability boundaries in the stability chart, Hopf bifurcations take place leading to self-excited vibrations in the corresponding nonlinear PDE-IDE system. The appearing periodic and quasi-periodic oscillations were found by numerical simulations and the stability chart was confirmed numerically. A constructed experimental rig allowed to detect the periodic and quasi-periodic oscillations experimentally and to determine an experimental stability chart The experimental stability chart confirmed the importance of the non-zero relaxation length. In previous experimental studies where zero relaxation length was considered [START_REF] Takács | Dynamics of Rolling of Elastic Wheels[END_REF]Takács and Stépán, to appear 2007), the results showed large deviation from theoretical predictions for large towing speed and long caster. The appearance of quasi-periodic self-excited oscillations at low towing speed and short caster also validates the delayed shimmy model. To explain these oscillations it is essential to describe appropriately the 'motion of the contact line'. Single degree-of-freedom models with creep force approximation cannot predict quasi-periodic behaviour because they omit the dynamics within the contact region. Nevertheless, whether a towing speed is considered to be large or small depends on the natural frequency of the system, and also on the length of the contact patch as it is expressed by the dimensionless towing speed.

Note that the proper bifurcation analysis of the nonlinear system has not been carried out yet. This can be very complicated in the studied infinite dimensional system in particular in the vicinity of the double Hopf bifurcation points. In order to resolve this problem one needs to consider the nonlinear terms and use either normal form calculations [START_REF] Campbell | Analytical and symbolically-assisted investigations of Hopf bifurcations in delay-differential equations[END_REF][START_REF] Orosz | Hopf bifurcation calculations in delayed systems[END_REF] or numerical continuation techniques [START_REF] Engelborghs | dde-biftool v. 2.00: A Matlab package for bifurcation analysis of delay differential equations[END_REF][START_REF] Szalai | Continuation of bifurcations in periodic delay-differential equations using characteristic matrices[END_REF]. For example, carrying out these calculations for nonlinear point contact models with rigid tyre, strong subcritical behaviour was found [START_REF] Stépán | Chaotic motion of wheels[END_REF]Takács et al., 2007). Subcriticality results in small-amplitude unstable oscillations around the stable stationary rolling and also predicts bistability between stationary rolling and large-amplitude oscillations.

Further development of our model is possible by involving sliding of the tyre at the rear part of the contact region. The appearing friction force dissipates energy and may allow us to have a better match between simulations and experiments (Takács and Stépán, to appear 2007). However, allowing the wheel to slide may also lead to very complicated (e.g., chaotic) large-amplitude vibrations in some parameter regimes [START_REF] Stépán | Chaotic motion of wheels[END_REF]Takács et al., 2007). We consider these problems for future research directions.

Fig. 2 .

 2 Fig. 2. Stretched string tyre model. The deformation of the central line of the tyre is shown in 3 dimension in panel (a) and projected to the ground, i.e., to the (x, y)-plane in panel (b). The relaxation length σ is identified corresponding to formula (1).

  33) satisfy (31) with the dimensionless angular frequency ω. Due to this bifurcation, self-excited vibrations may appear in the corresponding nonlinear system around the stationary rolling motion with dimensional frequency f = ωv/(4aπ) = ωV ω n /(2π) in Hertz. Consequently, travelling waves propagate backward along the contact line with dimensional wave length v/f = 4aπ/ω. The stability boundaries are determined in the parameter space by the substitution of the characteristic exponent λ 1 = iω into (31) and by the separation and imaginary parts: Re D(iω; μ cr ) = 0 , Im D(iω; μ cr ) = 0 . (34)

Fig. 3 .

 3 Fig. 3. Stability charts in the (V, L)-plane for the undamped system (ζ = 0) for different values of the dimensionless relaxation length Σ. Stable regions are shaded. The stability boundary curves refer to Hopf bifurcations, their intersections refer to co-dimension two Hopf bifurcations.

Fig. 4 .

 4 Fig. 4. Stability charts in the (V, L)-plane for different values of the damping ratio ζ with fixed dimensionless relaxation length Σ = 1.8. Stable regions are shaded. The stability boundary curves refer to Hopf bifurcations, their intersections refer to co-dimension two Hopf bifurcations. Thin curves correspond to the stability boundaries of the undamped system (ζ = 0) for the same Σ = 1.8 (see Fig. 3(d)). In panel (c) the crosses A-D show the parameter values used for the numerical simulations in Fig. 5.

Fig. 5 .

 5 Fig. 5. Numerical simulation results for the points marked by red crosses in Fig. 4(c).In each panel the time profile for the caster angle ψ(t) is shown and the spatial distribution of the lateral deformation q * (x) = q(x, t * ) is depicted at time t * . In panels (a) and (c) the stationary rolling motion is linearly stable, while in panels (b) and (d) this motion is unstable and stable oscillating (shimmy) motions can be observed.

Fig. 6 .

 6 Fig. 6. The experimental rig on the stiffened conveyor belt.Considering cases A and C the appearing dimensional frequency and wavelength are close to f = ωV ω n /(2π) = 0.93 f n and 4aπ/ω = 3.22 (2a), respectively, where the dimensionless angular frequency ω = 1.95 is obtained from the neighbouring Hopf curve. Similarly, in cases B and D the appearing dimensional angular frequencies are close tof 1 = ω 1 V ω n /(2π) = 0.27 f n and f 2 = ω 2 V ω n /(2π) = 1.03 f n ,and the appearing wavelengths are close 4aπ/ω 1 = 3.85 (2a) and 4aπ/ω 2 = 1.01 (2a), respectively. Here the dimensionless angular frequencies ω 1 = 1.63 and ω 2 = 6.20 belong to the neighbouring co-dimension two Hopf point.

Fig. 7 .

 7 Fig. 7. Measuring tyre parameters. Panel (a) shows the static lateral load on the framed wheel. Panel (b) shows the measurements of contact length 2a and relaxation length σ of the tyre contacting a transparent plastic plate.

Fig. 8 .

 8 Fig.8. Panel (a) compares the experimentally identified stability boundary (piecewise smooth line with circles) with the theoretical stability boundaries (smooth lines) on the (V, L)-plane. The experimentally identified stable region is shaded. For parameter values at the point D (marked by cross) quasi-periodic oscillations were found by numerical simulation and by experiment, too.

  Fig. 8 compares this measured stability chart with

  good match to the theoretical one.

Acknowledgments:

The research of D.T. and G.S. was supported by the Hungarian National Science Foundation under grant no. OTKA T043368. G.O. acknowledges the discussions with Andrew Gilbert on numerical methods for PDEs. The authors thank the reviewers for the detailed and useful comments.


S Wookey

H Kaneshima

T Inoue

P Yamashina

Harjadi

J O S Hammond
email: j.hammond@bristol.ac.uk

J Wookey

S Kaneshima

H Inoue

T Yamashina

P Harjadi

Systematic variation in anisotropy beneath the mantle wedge in the Java-Sumatra subduction system from shear-wave splitting

Systematic variation in anisotropy beneath the mantle

wedge in the Java-Sumatra subduction system from shear-wave splitting

shear-wave splitting in local and teleseismic data from 12 broadband stations across Sumatra and Java to study the anisotropic characteristics of this subduction system, which can provide important constraints on dynamical processes involved. Splitting in S-waves from local earthquakes between 75-300 km deep show roughly trench parallel fast directions, and with timelags 0.1-1.3 s (92% ≤0.6 s). Splitting from deeper local events and SKS, however, shows larger time-lags (0.8-2.0 s) and significant variation in fast direction. In order to infer patterns of deformation in the slab we apply a hybrid modelling scheme. We raytrace through an isotropic subduction zone velocity model, obtaining event to station raypaths in the upper man-A c c e p t e d M a n u s c r i p t tle. We then apply appropriately rotated olivine elastic constants to various parts of the subduction zone, and predict the shear-wave splitting accrued along the raypath. Finally, we perform grid searches for orientation of deformation, and attempt to minimise the misfit between predicted and observed shear-wave splitting. Splitting from the shallow local events is best explained by anisotropy confined to a 40 km over-riding plate with horizontal, trench parallel deformation. However, in order to explain the larger lag times from SKS and deeper events, we must consider an additional region of seismic anisotropy in or around the slab. The slab geometry in the model is constrained by seismicity and regional tomography models, and many SKS raypaths travel large distances within the slab. Models placing anisotropy in the slab produce smaller misfits than outside for most stations. There is a strong indication that inferred flow directions are different for sub-Sumatran stations than for sub-Javanese, with >60 • change over ∼375 km. The former appear aligned with the subduction plate motion, whereas the latter are closer to perpendicular, parallel to the trench direction. There are significant differences between the slab being subducted beneath Sumatra, and that beneath Java: age of seafloor, maximum depth of seismicity, relative strength of the bulk sound and shear-wave velocity anomaly and location of volcanic front all vary along the trench. We speculate, therefore, that the anisotropy may be a fossilised signature rather than due to contemporary dynamics.

Key words: Seismic anisotropy; Subduction zones; Indonesia; S-wave At the Java-Sunda trench the Indian-Australian plates are being subducted beneath the Eurasian plate. This is one of the most seismically active regions in the world, with earthquakes occurring from the surface to the base of the transition zone. This makes the region an ideal natural laboratory to study subduction dynamics using seismic imaging techniques. Here, we present a study of seismic anisotropy (the variation of seismic wavespeed with direction) beneath the Java-Sunda subduction system, providing some insights into the deformation history and active processes in a subduction setting.

As a seismic shear-wave propagates through an anisotropic medium it is split into two quasi shear-waves, one with particle motion in the direction of fast seismic wavespeed, and another perpendicular to this. This results in the two shear-waves travelling at different speeds and arriving at a seismic station separately; a phenomenon called shear-wave splitting (see [START_REF] Crampin | Shear-wave polarizations near the North Anatolian Fault, II, Interpretation in terms of crack-induced anisotropy[END_REF][START_REF] Silver | Shear-wave splitting and subcontinental mantle deformation[END_REF]. At a 3-component seismic station it is possible to measure these separate shear-waves and determine the fast direction of wave polarisation and the time between the shear-waves (see [START_REF] Silver | Seismic anisotropy beneath the continents: Probing the depths of geology[END_REF][START_REF] Savage | Seismic anisotropy and mantle deformation: What have we learned from shear wave splitting?[END_REF]Kendall, 2000, for reviews).

In many cases anisotropy in the upper mantle is assumed to be caused by the lattice preferred orientation (LPO) of olivine. Typically, the olivine fast axis (a-axis) aligns in the direction of upper-mantle flow (e.g., [START_REF] Babuska | Seismic Anisotropy in the Earth[END_REF][START_REF] Mainprice | The seismic anisotropy of the Earth's mantle: from single crystal to polycrystal[END_REF][START_REF] Mainprice | Seismic anisotropy of the deep Earth from a mineral and rock physics perspective[END_REF]. As a result shearwave splitting can provide direct information about dynamic processes, such as mantle flow, and accumulated strain due to previous deformation events A c c e p t e d M a n u s c r i p t which have 'frozen' a source of anisotropy into the lithosphere beneath a seismic station.

This technique has been applied to many subduction zone settings around the world, both using local slab related events, and teleseismic events (e.g., [START_REF] Fouch | Mantle anisotropy beneath northwest Pacific subduction zones[END_REF][START_REF] Long | The subduction zone flow field from seismic anisotropy: A global view[END_REF]. The results are far from simple with fast directions appearing trench parallel, perpendicular and oblique, often at the same subduction zone [e.g., New Zealand, [START_REF] Morley | Mantle wedge anisotropy in the Hikurangi subduction zone, central North Island, New Zealand[END_REF], Japan [START_REF] Nakajima | Shear-wave polarization anisotropy and A c c e p t e d M a n u s c r i p t subduction-induced flow in the mantle wedge of northeastern Japan[END_REF], South America [START_REF] Helffrich | A teleseismic shear-wave splitting study to investigate mantle flow around South America and implications for plate-driving forces[END_REF] and the Carribean (Pinero-Felicangeli and [START_REF] Pinero-Felicangeli | Sub-slab mantle flow parallel to the Caribbean plate boundaries: Inferences from SKS splitting[END_REF]]. This suggests a complex region of anisotropy and many different mechanisms have been invoked to explain these results.

Deformation in the overriding plate and return flow in the mantle wedge are discussed as causing trench parallel and perpendicular fast directions for New Zealand [START_REF] Morley | Mantle wedge anisotropy in the Hikurangi subduction zone, central North Island, New Zealand[END_REF]. The presence of water has been suggested to change the alignment mechanism of olivine, thus producing fast directions perpendicular to the direction of maximum stress [START_REF] Jung | Water-induced fabric transitions in olivine[END_REF] and this has been invoked to explain trench parallel splitting beneath Japan [START_REF] Nakajima | Shear-wave polarization anisotropy and A c c e p t e d M a n u s c r i p t subduction-induced flow in the mantle wedge of northeastern Japan[END_REF], although the authors also acknowledge that trench parallel flow, caused by along strike dip variations might also explain these results. Evidence for deeper sources of anisotropy has also been described to explain splitting seen in teleseismic shear-wave arrivals, for example slab anisotropy in New Zealand [START_REF] Brisbourne | Anisotropic structure of the Hikurangi subduction zone, New Zealand-integrated interpretation of surface-wave andbody-wave observations[END_REF], sub-slab anisotropy due to toroidal flow beneath the slab in the Carribean [START_REF] Pinero-Felicangeli | Sub-slab mantle flow parallel to the Caribbean plate boundaries: Inferences from SKS splitting[END_REF] and basal drag beneath the slab in South America [START_REF] Helffrich | A teleseismic shear-wave splitting study to investigate mantle flow around South America and implications for plate-driving forces[END_REF].

Page 5 of 47 Another model invoked to explain global observations of subduction zone anisotropy suggests preferentially-oriented hydrated faults in the subducting slab [START_REF] Faccenda | Fault-induced seismic anisotropy by hydration in subducting oceanic plates[END_REF]. During the subduction process, faulting in the downgoing plate allows water to migrate into the crust and mantle. This results in the formation of sheets of silicates in the top few kms of the slab, inducing a shape preferred orientation, resulting in high degrees of anisotropy.

Superimposed on this, an alignment of serpentinized rocks will cause an additional crystallographic preferred orientation, thus giving the potential for large degrees of shear-wave splitting from a thin layer at the top of the slab [START_REF] Faccenda | Fault-induced seismic anisotropy by hydration in subducting oceanic plates[END_REF].

The variability in shear-wave splitting results, and the wide variety of A c c e p t e d M a n u s c r i p t mechanisms described to explain them shows that subduction zone settings are very complicated. Possible regions of anisotropy exist in the crust, lithosphere and mantle wedge, and in or below the slab. In this study we try to understand the problem in more detail by forward modelling our shear-wave splitting results with realistic subduction zone models.

Data and methodology

We use data from 9 stations of the JISNET (Japan Indonesia Seismic Network) array [START_REF] Ohtaki | Sunarjo, Broadband Seismic Network in Indonesia-JISNET[END_REF] to determine the anisotropic characteristics beneath Sumatra and Java. The stations were deployed in 1998 and some stations are still recording. We use data up until the end of 2003. We also use 2 stations from the GEOFON array (UGM, 9 years, GSI, 3 years) and 1 from the Ocean hemisphere network project (PSI, 13 years), a station used by [START_REF] Long | The subduction zone flow field from seismic anisotropy: A global view[END_REF] to investigate subduction zone shear-wave splitting. Station locations are shown in Figure 1.

We employ earthquakes with distances between 85 • -140 • , events ideal for SKS/SKKS-wave splitting analysis. A total of 47 earthquake-station paths provided SKS/SKKS-phases of sufficient quality to enable splitting to be estimated (2% of data used). All SKS/SKKS data were filtered between 0.05-0.3 Hz. We also measure splitting in S-waves from local earthquakes where the incidence angle is less than 45 • [i.e., within the shearwave window [START_REF] Evans | Effects of the free surface on shear wavetrains[END_REF]]. Earthquake locations from the catalogue of [START_REF] Engdahl | Global teleseismic earthquake relocation with improved travel times and procedures for depth determination[END_REF] S1-S4. All local earthquake data were filtered between 0.1-1 Hz. The overlap between SKS/SKKS-and local S-wave filter bands helps to minimise frequency dependent effects between the two data sets.

We estimate shear-wave splitting in both the core phases and S-phases from local earthquakes using the method of Teanby et al. (2004), which is based on the methodology of [START_REF] Silver | Shear-wave splitting and subcontinental mantle deformation[END_REF]. It uses a grid search over fast direction (φ) and time-lag (δt), minimising the second eigenvalue of the covariance matrix for the particle motion around a given time window. The advantage of the Teanby et al. ( 2004) method is that it automates window selection, making 100 measurements around the relevant phase, and applying cluster analysis to identify the most stable result. Figure 2 displays a typical result of a local S-and teleseismic SKS-wave splitting result at PSI.

Results

From our initial analysis we reject results with 1σ errors greater than 0.4 s in δt and 30 • in φ. Furthermore, the polarisation of the shear-wave after the anisotropy correction is applied (source polarisation, hereafter spol) is compared with the back azimuth for SKS arrivals. For radially polarised phases (i.e a P-S conversion at the core/mantle boundary) these should be similar, thus we reject results with differences >30 • . Local earthquakes are not radially polarised as they have not been converted to a P-wave at any point along their ray path. However some earthquakes are large enough to have a previously determined focal mechanism (we use the Global CMT and USGS catalogue), from which the initial polarisation can be calculated.

A c c e p t e d M a n u s c r i p t

For these earthquakes the spol measurement and initial polarisation are compared to provide an extra quality assessment (10% of local splitting measurements included this check). Again spol and the initial polarisation must agree within 30 • to be included. Those results coming from earthquakes without an available focal mechanism are assumed to be correct. Results are shown in Figure 3 and Tables S1-S4. We discuss the results for Sumatra and Java separately in the following sections. [START_REF] Gudmundsson | A regionalized upper mantle(RUM) seismic model[END_REF]). As a result splitting is estimated in shear-waves generated from earthquakes in the ∼100 km-200 km depth range. One station, GSI, is located closer to the trench, and thus provides estimates of splitting from earthquakes in the depth range ∼30-80 km. All these splitting results show little variation with depth (Figure 4). Time-lags for all stations, including GSI, are generally between 0.2-0.4 s, and fast directions roughly trench, parallel, with some scatter, and parallel to the strike of the Sumatran fault (Figures 3 and4). This indicates that the source of this anisotropy is in the over-riding plate and that the mantle wedge above the slab is largely isotropic.

SKS/SKKS-wave splitting

The stations located above the 150 km slab contour (BSI, PSI, PPI, KSI and KOTA) show ∼1 s of splitting, with a general north-south trend. How-A c c e p t e d M a n u s c r i p t ever some variation in the fast direction is evident, particularly at PSI (Figure 3). The station close to the trench, GSI, is markedly different with a large time-lag of 1.8 s and a trench perpendicular fast direction of 52.5 • . The obvious discrepancy between the time-lags for local and teleseismic splitting suggests that a deeper source of anisotropy, within or beneath the slab, must be present.

Java

Local S-wave splitting

Local earthquakes, within the shear-wave window, have been recorded down to depths of 700 km at most stations in Java. Splitting results from these events display considerable variation in both fast direction and time lag on individual stations as well as from different stations (Figure 3). Most stations show both trench parallel and trench perpendicular fast directions. The trench perpendicular results are in general from intermediate depth earthquakes (250-350 km), with both shallower and deeper earthquakes having more trench parallel fast directions. Some deep events produce large delay times and fast directions oblique to the strike of the trench (UGM, BRB).

These seem isolated to one region of the subduction zone, with adjacent splitting measurements from nearby earthquakes showing much smaller delay times, albeit with similar fast directions. Splitting at BMI is very complicated, possibly because it lies above a part of the slab which bends beneath Sulawesi.

A c c e p t e d M a n u s c r i p t

SKS/SKKS-wave splitting

In general SKS-wave results have an east-west fast direction, with ∼1.6 s of splitting (Figure 3). There does seem to be some subtle variation in the splitting results with the estimated fast directions occurring in two groups, ∼75 • and ∼105 • . The magnitude of splitting is again larger than the splitting seen in local events, except for a few deep earthquakes at UGM, and BRB.

Again, similarly to the results at Sumatra, this suggests a slab or sub-slab source of anisotropy.

Modelling

It is evident that more than one region of anisotropy exists beneath Sumatra and Java. Shallow earthquakes (<300 km deep) beneath Sumatra and Java show little variation in time-lag with depth, and show consistent trench parallel fast directions. This suggests that the mantle wedge is isotropic, and that the splitting observed from these earthquakes is accrued in the overriding plate. However larger time-lags observed in SKS/SKKS splitting at Sumatra and Java, and splitting from deeper (>300 km deep) earthquakes beneath Java suggest that anisotropy must exist in the slab, or sub-slab as well.

To understand the anisotropy in more detail we have developed a technique based on ray-tracing through an inhomogeneous subduction zone model [START_REF] Guest | Modelling seismic waveforms in anisotropic inhomogeneous media using ray and maslov asymptotic theory: Applications to exploration seismology[END_REF]. We build a model consisting of four regions: layer 1, the over-riding plate and mantle wedge; layer 2, a thin 10 km layer above the slab; layer 3, the subducting slab; and layer 4, the sub-slab mantle (Figure 5). The surface of the slab is calculated Page 11 of 47 A c c e p t e d M a n u s c r i p t using the regionalized upper mantle (RUM) seismic model slab contours [START_REF] Gudmundsson | A regionalized upper mantle(RUM) seismic model[END_REF] which is consistent with a relocated local seismicity study [START_REF] Fauzi | Lateral variation in slab orientation beneath Toba Caldera, northern Sumatra[END_REF] and deep events from the USGS earthquake catalogue. The model is also consistent with regional tomography models for Indonesia [START_REF] Gorbatov | Joint bulk-sound and shear tomography for Western Pacific subduction zones[END_REF]. Uncertainty in these slab models could introduce error into these results, however recently published slab models of Syracuse and Abers (2006) show only minor changes (∼20 km shift east) in Sumatra-Java subduction zone to those calculated by [START_REF] Gudmundsson | A regionalized upper mantle(RUM) seismic model[END_REF].

We use ATRAK, a ray-tracer capable of tracking seismic rays in multilayered 3-D media with general anisotropy [START_REF] Guest | Modelling seismic waveforms in anisotropic inhomogeneous media using ray and maslov asymptotic theory: Applications to exploration seismology[END_REF]. This technique has been utilised to investigate anisotropy in other subduction zones [e.g., Kurils [START_REF] Kendall | Seismic modelling of subduction zones with inhomogeneity and anisotropy-I. Teleseismic P-wavefront tracking[END_REF] and New Zealand [START_REF] Brisbourne | Anisotropic structure of the Hikurangi subduction zone, New Zealand-integrated interpretation of surface-wave andbody-wave observations[END_REF]]. Our model differs from previous studies by the fact that we trace rays through an isotropic model, assuming that anisotropy only minimally distorts the raypath (a reasonable assumption for relatively weakly anisotropic media). We trace rays for each station separately (except UGM and SWH which are very close together, and so are incorporated into one model). We appeal to reciprocity and shoot rays from the station, saturating the model (>12,000 rays per station, Figure 5). We then match each event to its nearest ray, thus defining the raypath taken from the earthquake to the station. For local earthquakes we use the latitude, longitude and depth to match the ray, and for SKS/SKKS phases we calculate the piercing point of the phase at a depth of 800 km beneath the station [through the A c c e p t e d M a n u s c r i p t AK135 1-D velocity model [START_REF] Kennett | Constraints on seismic velocities in the Earth from traveltimes[END_REF]], and use this information as the event location (Figure 5). This ray-path for each station-event pair is then broken down into a length, time and angle spent in each part of the model. With this information we can then impose an anisotropy on a specified region of the model (layers 1-4), represented by an appropriately rotated elastic tensor (we use olivine elastic constants from [START_REF] Abramson | The elastic constants of San Carlos olivine to 17 GPa[END_REF]). Anisotropy is confined to the olivine stability field (i.e. <410 km), and can be restricted to a particular depth range within a layer.

For each contiguous section of the raypath we can use the Christoffel equation ( Červenỳ, 2001) to calculate the fast and slow shear-wave velocities and particle motion, thus, along with the distance in the layer, we can calculate a time-lag and fast direction through the layer.

We then perform a grid search over orientation (γ, ∆γ=5 • ), and anisotropic strength (by Voigt-Reuss-Hill dilution of the elastic tensor; denoted volume fraction aligned -VFA, ∆VFA=0.01-0.05).

The angle γ is the orientation of the olivine a-axis in a plane which is either horizontal (shallow models) or inclined at the average dip of the subducting slab (deep models). The b-axis is held normal to this plane. For each node in the grid we calculate δt and φ for each event-station raypath in the data. The model is then evaluated as the summed misfit between these and the real measurements (for δt and φ separately, and for a combination). Typically, VFA is poorly constrained, so is not interpreted in great detail, instead we primarily analyse the best fitting orientation γ, illustrated by the density rotorgram in Figure 6.

Page 13 of 47 This can be performed for each layer in the model. Four models are run in total, the first constraining the layer 1 anisotropy, and the last three incorporate this result with layer 2, 3 and 4 respectively. Multiple layers of anisotropy can be included, though only one may vary, as a result we can not formely asses trade-off between these layers. Where two layers of anisotropy are included we combine the splitting operators using the n-layer splitting equations of [START_REF] Silver | The interpretation of shear wave splitting[END_REF], and the initial source polarisation from the splitting measurements. An example result is shown in Figure 6 for station GSI with anisotropy in layer 1 only.

The locations of stations and earthquakes used in this study are ideal to image the various parts of the subduction system. We image paths in the sub-slab, slab, supra-slab, mantle wedge and over-riding plate. This enables us to go one step further than typical shear-wave splitting studies and place constraints on the location of the anisotropy, and thus better understand the dynamics of Java-Sunda subduction.

The modelling scheme outlined here allows us to extract information on which part of the subduction system each event is sampling. Therefore, we are able to investigate further the observation made in section 3, that the mantle wedge appears largely isotropic. Figure 7 shows the variation in delay time as a function of path length through various parts of the model. It is evident that there is no clear relationship between splitting delay times and path length in the mantle wedge. In fact, those data with paths almost exclusively in the wedge display the smallest delay times (<0.4 s, except two SKS events at PSI which, in our model, never sample the region below the A c c e p t e d M a n u s c r i p t wedge). An estimate can be placed on the amount of anisotropy in a 40 km thick over-riding plate and the mantle wedge below. Figure 7a shows a linear fit to the data with paths exclusively in the wedge (removing the two SKSwave arrivals at PSI). The line has equation δt = 0.0007l + 0.1306, where l is the path length in the mantle wedge and over-riding plate. Assuming an average S-wave velocity of 4.5 km/s, this equates to an anisotropy of 1.8% in a 40 km thick plate, and 0.3% anisotropy in the mantle wedge below this plate. This shows that the anisotropy in the over-riding plate is dominating the shear-wave splitting. In comparison, some correlation is evident between splitting delay time and path length spent in the region beneath the mantle wedge, but above 410 km (layers 2-4 in our model) (Figure 7b). The splitting results can be explained by a model which has 20% aligned olivine crystals along this path length. This observation supports the idea that the mantle wedge is largely isotropic, and that anisotropy is largely confined to the over-riding plate, slab and sub-slab region. To further constrain the location of the anisotropy we split the modelling into two steps, shallow and deep anisotropy. To constrain the characteristics of the shallow anisotropy we use events with depths <300 km, and to model deeper anisotropy we use events with depths >300 km and SKS-events.

Shallow anisotropy

We impose anisotropy in the uppermost 40 km of layer 1 of the model, simulating the overriding plate lithosphere. The choice of this thickness is fairly arbitrary (i.e. a layer half as thick with twice the anisotropy strength will fit the data just as well), but we base our model on the limit Page 15 of 47 A c c e p t e d M a n u s c r i p t of our splitting results (the shallowest events at GSI are ∼30 km). In these cases the olivine a-axis is horizontal and the misfit is calculated only using events shallower than 300 km.

Table 1 and Figure 8 show that the best fitting anisotropy orientations at all stations are well constrained and show trench parallel orientations. BRB is not modelled as we have too few splitting measurements. Errors in depth locations could have an effect on our estimated VFA, however we expect this effect to be minimal due to the fact that we impose a 40 km thick layer of anisotropy at the top of the model. Events which are mislocated from below 40 km to above 40 km, and events which are mislocated to below 40 km from above 40 km will cause the VFA to be over and underestimated respectively. This is only a problem at GSI, where event depths range from 30-80 km. At other stations only one event occurs shallower than 70 km (SWH, 33 km).

Deep anisotropy

To fit local slab events with depths >300 km, and SKS/SKKS phases we test 3 classes of model where anisotropy is confined to a thin layer above the slab (layer 2), the slab (layer 3) or the sub-slab (layer 4). We impose an olivine anisotropy with its a-axis oriented (γ) in a plane inclined at the average dip of the subduction zone (baxis is normal to the plane). In all deep models we also impose the best fitting shallow anisotropy for each station (determined above). Thus, each model has two anisotropic regions, one fixed in orientation and strength (shallow; layer 1) and one varying (deep;

A c c e p t e d M a n u s c r i p t layers 2-4). Figures 8 and9 andTable 2, and it is evident that it is hard to confidently discriminate between the slab and sub-slab models (layer 3 and 4), with both giving compatible anisotropy orientations.

Results are shown in

Anisotropy in a thin layer at the top of the slab (layer 2) consistently has the largest misfit, and requires a very high degree of alignment, suggesting this is a less likely cause of our observations (Figure 9 and Table 2). For both Java and Sumatra the misfit is considerably lower for a model where anisotropy is confined to the slab (assuming a thickness of 100 km). A robust feature in all models (with the exception of GSI) is that the deeper anisotropy orientation beneath Sumatra is more north-south, and rapidly changes at the Sunda strait becoming more east-west beneath Java (Figures 8 and9). GSI shows a trench perpendicular anisotropy orientation, suggesting that a different mechanism is causing anisotropy at the trench compared to that beneath the volcanic/back-arc region (Figure 8). GSI has the longest path in the sub-slab region. The anomolous trench perpendicular anisotropy orientation, and large time delays at this station may suggest that sub-slab anisotropy is dominant in the forearc. BRB and KHK do not have enough observations to model. Mislocated events for deep and SKS/SKKS-wave events will have minimal effects on their raypaths.

Discussion

From the modelling it is evident that anisotropy in the over-riding plate is required to explain the splitting observations from events <300 km deep at Java and Sumatra. The modelled anisotropy orientation is near paral-Page 17 of 47 A c c e p t e d M a n u s c r i p t lel to the strike of the trench. The simplest explanation for this anisotropy is vertically aligned cracks associated with deformation in the over-riding Eurasian plate, with a largely isotropic mantle wedge below. This observation is supported by a recent local earthquake tomography study for central Java, where high P-wave anisotropy (7-10%) is observed in the crust and lithospheric mantle [START_REF] Koulakov | Anisotropic structure beneath central Java from local earthquake tomography[END_REF]. Similar observations of overriding plate anisotropy and isotropic mantle wedges have been seen at New Zealand [START_REF] Morley | Mantle wedge anisotropy in the Hikurangi subduction zone, central North Island, New Zealand[END_REF], the Carribean (Pinero-Felicangeli and [START_REF] Pinero-Felicangeli | Sub-slab mantle flow parallel to the Caribbean plate boundaries: Inferences from SKS splitting[END_REF] and South America [START_REF] Polet | Shear wave anisotropy beneath the Andes from the BANJO, SEDA and PISCO experiments[END_REF]. The fact that the wedge is isotropic suggests that either 2D corner flow caused by viscous coupling of the slab and wedge is of similar magnitude as lateral flow from trench migration, or that they are both weak in this region, thus meaning no coherent LPO can be developed [START_REF] Long | The subduction zone flow field from seismic anisotropy: A global view[END_REF]. This observation is in drastic contrast to other subduction zones such as Japan, where Hiramatsu et al. [START_REF] Long | The subduction zone flow field from seismic anisotropy: A global view[END_REF] suggest that 2D corner flow is causing a large anisotropy. In all these regions they observe large δt ∼1-2 s, much larger than we observe beneath Java and Sumatra. Models of deeper anisotropy, constrained by deep subduction events, and SKS/SKKS show that a slab source of anisotropy is most likely, although a sub-slab source can not be ruled out and gives similar best fitting anisotropy orientations. However, in all these models a sharp change in anisotropic characteristics is observed between Sumatra and Java (Figure 8). This change A c c e p t e d M a n u s c r i p t is also evident in the raw data (Figure 10). SKS fast direction estimates for the same event recorded at KOTA, the southernmost Sumatran station, and LEM the westernmost Javan station vary by ∼60 • over ∼375km (Figure 10). This change in anisotropic characteristics is mirrored by changes in other geophysical parameters close to the Sunda Strait. The most obvious change near this location is in the age of the downgoing slab (Figure 1). Magnetic anomalies show a change in age from ∼60Ma to ∼100Ma from southern Sumatra to Western Java [START_REF] Sdrolias | Controls on back-arc basin formation[END_REF]. The sharpness, and exact location of this boundary is unclear, especially at depth in the subducted slab. Syracuse and Abers (2006) observed that this age boundary co-incides with a change in volcanic characteristics. They found that the volcanic front in Sumatra is located above the 90 km slab contour, with an abrupt change close to the Sunda Strait where volcanism is present above the 150 km slab contour. This change occurs over <150 km, similarly to the abrupt change in anisotropy, and also correlates with the abrupt reduction in deep (>500 km) seismicity beneath Sumatra (e.g., [START_REF] Engdahl | Global teleseismic earthquake relocation with improved travel times and procedures for depth determination[END_REF], although no change in the shape of the slab is observed. Syracuse and Abers (2006) also correlate this change in the location of the volcanic front with variation in the ratio of potassium to silica (Wheller et al., 1987), but the reasons for this are unclear (Syracuse and Abers, 2006). Regional tomography profiles of [START_REF] Gorbatov | Joint bulk-sound and shear tomography for Western Pacific subduction zones[END_REF] The alignment of the deep anisotropy orientations at all stations beneath Sumatra, except GSI, roughly matches with the APM of the subducting slab and oblique to the orientation of the trench. This would suggest that in this case shear associated with the subducting plate moving through the mantle is aligning sub-slab olivine, and causing APM aligned anisotropy.

However, beneath Java, the opposite is true. The deep anisotropy orientation is perpendicular to the APM of the subducting slab. The slab between Sumatra and Java appears continuous (Syracuse and Abers, 2006), so it seems unlikely that basal drag would create an anisotropy beneath southern Sumatra, and not western Java <375 km away. Thus, basal drag does not seem a likely cause of anisotropy. [START_REF] Long | The subduction zone flow field from seismic anisotropy: A global view[END_REF] suggest that lateral flow beneath the slab, produced by trench migration, can induce a trench parallel anisotropy orientation, and this is observed beneath Java. However, this would predict trench parallel anisotropy orientations beneath Sumatra, and we observe trench oblique orientations. Again, this change from lateral flow beneath western Java, to an oblique flow beneath southern Sumatra seems unlikely over such a short distance.

A c c e p t e d M a n u s c r i p t

Water in the mantle

Jung and Karato ( 2001) have suggested that the presence of water can make the b-axis of olivine align with the flow direction, manifesting in a anisotropy orientation perpendicular to the direction of flow. Thus a change in the amount of water beneath Java and Sumatra could induce a change from APM parallel anisotropy orientations to APM normal anisotropy orientations. The most likely region of water rich mantle is above the slab, in the mantle wedge. We observe that the mantle wedge in this region is largely isotropic, with little variation in time-lags with depth, suggesting that this is an unlikely scenario to cause the sharp change in anisotropy.

Fossil anisotropy

Figure 5 shows that the deep and teleseismic events travel some distance in the subducting slab, and thus anisotropy in the slab will manifest itself in the observed splitting results. Interestingly, the sharp change in anisotropic characteristics between Java and Sumatra co-incides with major changes in slab characteristics.

The age of the subducting plate beneath Java and Sumatra varies between 100 Ma in the north to 43 Ma at 2 • S and up to 160 Ma in the east [START_REF] Sdrolias | Controls on back-arc basin formation[END_REF]. This transition in ages occurs near the region where we see an abrupt change in anisotropy with the plate subducting beneath Sumatra and Abers, 2006). The abrupt changes in age and seismic properties suggests that differences between the Javan and Sumatran slabs may be the cause for the abrupt change in anisotropic characteristics. [START_REF] Faccenda | Fault-induced seismic anisotropy by hydration in subducting oceanic plates[END_REF] show that preferentially oriented hydrated faults in the uppermost part of the slab can produce a high degree of shear-wave splitting. A change in age between the slab beneath western Java and southern Sumatra may effect the efficacy of fracture generation in the slab, and thus hydrated faults may change in style or quantity across this region. However, our models indicate that a thin highly anisotropic layer at the top of the slab poorly fits our results, mainly due to the fact that rays spend such a small amount of time in this region. If this layer behaves as a strong waveguide, a feature not modelled by our simple subduction models, then this may explain the observed shear-wave splitting.

Another possible location for a fossil anisotropy is in the slab lithosphere.

As oceanic lithosphere is formed and moves over the asthenosphere it creates a large anisotropic signature (e.g. at the East Pacific Rise Wolfe and Solomon, 1998;[START_REF] Harmon | Variations in shearwave splitting in young Pacific seafloor[END_REF], in the direction of plate motion. These anisotropic characteristics can be frozen into the lithosphere, and may remain present in the subducted slab beneath Java and Sumatra. This mechanism has been suggested to be responsible for at least part of the observed shear-wave splitting seen at the Marianas [START_REF] Fouch | Shear wave anisotropy in the Mariana subduction zone[END_REF]. Also, [START_REF] Brisbourne | Anisotropic structure of the Hikurangi subduction zone, New Zealand-integrated interpretation of surface-wave andbody-wave observations[END_REF] show that significant anisotropy, with trench parallel anisotropy orientations, exists in the Hikurangi subduction zone, New Zealand. The Pacific plate in this region has a trench parallel age progres-A c c e p t e d M a n u s c r i p t sion [i.e. the age of the Pacific oceanic lithosphere changes in a direction parallel to the trench [START_REF] Sdrolias | Controls on back-arc basin formation[END_REF]]. This is similar to that seen beneath Java, suggesting that a fossilised anisotropic signature, frozen in as the plate formed, could explain anisotropy in this region. The oceanic lithosphere subducting beneath Java is 80-140 Ma, and Lithgow-Bertelloni and Richards (1998) show that plate motions for the Australian plate between 74-119Ma are generally east-west. The oceanic lithosphere subducting beneath Sumatra is 40-80 Ma, and Lithgow-Bertelloni and Richards (1998) show that the plate motions for the Indian plate at this time were northsouth. The shear-wave splitting change observed between Java and Sumatra could, therefore, be an effect of sampling these differing frozen anisotropic signatures.

Conclusion

Shear-wave splitting has been measured in SKS-wave arrivals and S-waves from local subduction zone earthquakes beneath the Java-Sunda subduction zone. Results show that the mantle wedge is mainly isotropic, with anisotropy largely confined to the over-riding plate. This has been further emphasised with modelling, which shows a best-fitting trench-parallel anisotropy orientation, suggesting anisotropy seen in splitting from local earthquakes is derived from anisotropy caused by vertically aligned cracks associated with deformation in the over-riding Eurasian plate. Splitting from deeper earthquakes, and SKS-wave arrivals have much larger delay times, suggesting a deeper source of anisotropy is present.

Modelling constrains the location of this anisotropy, with the most likely

A c c e p t e d M a n u s c r i p t anisotropic region being the slab. The best fitting anisotropy orientations in the slab were predominantly north-south beneath Sumatra (except GSI which lies above the trench, and shows trench perpendicular anisotropy orientationss, suggesting a different mechanism to that seen at the volcanic arc/back arc). This change in orientation occurs over less than 375 km and at the same location as changes in other geophysical characteristics such as slab age [START_REF] Sdrolias | Controls on back-arc basin formation[END_REF], bulk sound speed [START_REF] Gorbatov | Joint bulk-sound and shear tomography for Western Pacific subduction zones[END_REF] and location of the volcanic front (Syracuse and Abers, 2006), suggesting that a change in the slab characteristics is the reason for the change in anisotropy.

We speculate that the anisotropy has been frozen into the subducting lithosphere as it was formed, or moved over the asthenosphere. The subducting lithosphere beneath Java is older, and would have formed with eastward plate motions [START_REF] Lithgow-Bertelloni | The dynamics of Cenozoic and Mesozoic plate motions[END_REF], resulting in east-west anisotropy orientations, compared with the younger lithosphere subducting beneath Sumatra which formed when the Indian-Australian plate was moving north [START_REF] Lithgow-Bertelloni | The dynamics of Cenozoic and Mesozoic plate motions[END_REF], resulting in north-south anisotropy orientations.

Acknowledgments

We would like to thank Mike Kendall and George Helffrich for valuable discussion and help improving this paper. Also, we would like to thank Alexei Gorbatov for the use of his tomographic models. All the data used in this paper was provided by the Japan Indonesian Seismic Network (JISNET), IRIS, GEOFON and the Ocean Hemisphere Network Project, courtesy of [START_REF] Gripp | Current plate velocities relative to the hotspots incorporating the Nuvel-1 global plate motion model[END_REF], black arrows), slab contours at 100 km depth intervals [START_REF] Gudmundsson | A regionalized upper mantle(RUM) seismic model[END_REF], the Great Sumatra Fault, and Isochrons ( [START_REF] Müller | Digital isochrons of the world's ocean floor[END_REF]). 5) and b) path length above 410 km but beneath the mantle wedge (layers 2-4 in Figure 5) for all local and SKS splitting results. Red crosses in a) show splitting estimates for events which spend most of their path in the mantle wedge and over-riding plate (i.e. less than 10 km outside of this region). The red line shows the best fitting slope to the red crosses estimating an anisotropy of 1.8% in a 40 km thick over-riding plate, and 0.3% anisotropy in the mantle wedge below. The grey shaded region shows the range of potential splitting for an olivine fabric with 20% alignment (based on the minimum and maximum shear-wave anisotropy of the olivine elastic tensor). Note the poor correlation between the shaded area and the splitting results in a), and the consistently low delay times for events which spend most of their time in the mantle wedge and over-riding plate. The splitting results are explained much better by putting anisotropy in the slab or subslab region (panel b). This shows that the mantle wedge is largely isotropic. A) Best fitting orientations obtained from forward modelling. A) Olivine orientations best fitting the splitting results obtained from local events <300 km deep; these are predominantly trench parallel. B) Olivine orientations best fitting splitting results from local events >300 km deep and SKS/SKKS-wave splitting results. Note the rotation from northsouth orientations beneath Sumatra (GSI is an exception), and the east-west anisotropy orientations beneath Java. The complicated pattern at BMI reflects the scatter observed in the data. We speculate that this is due to complications arising from the slab bending beneath this station. 8 and Table 1) are included using the n-layer splitting equations of [START_REF] Silver | The interpretation of shear wave splitting[END_REF].

  Some attempts have been made to explain subduction zone shear-wave splitting observations with uniform global models.[START_REF] Long | The subduction zone flow field from seismic anisotropy: A global view[END_REF] compiled global observations of shear-wave splitting at subduction zones. They postulated that the largely trench parallel splitting observations which arise from sub-mantle wedge regions (i.e. slab or sub-slab) are dominated by lateral flow in the sub-slab region. This is induced by return flow due to trench migration.[START_REF] Long | The subduction zone flow field from seismic anisotropy: A global view[END_REF] explain mantle wedge anisotropy (above the slab) as being due to the interaction of lateral flow due to trench migration and corner flow driven by viscous coupling between the slab and the wedge. Such corner flow would induce a trench perpendicular fast direction in the splitting. The authors argue that the observed supra-slab trench parallel or trench perpendicular splitting results observed globally are due to the dominance of one or the other of these mechanisms.

  are used where possible, others are taken from the USGS catalogue, which could introduce some significant mislocation errors into our results. This results in 127 measurements of splitting Page 7 of 47 A c c e p t e d M a n u s c r i p t from local earthquakes (49% of data used). Individual results are shown in Tables

  1. Local S-wave splitting Seismic stations BSI, PSI, PPI, KSI and KOTA are located approximately above the 150 km contour of the subducting Indian plate (based on the slab contours of

  A c c e p t e d M a n u s c r i p t

(

  1998) observe two 75-175 km deep anisotropic zones that are caused by melt filled cracks, Tonga where Smith et al. (2001) suggest that along strike flow in the mantle related to slab rollback and the Samoan plume is inducing an anisotropy and the Aleutians where

  show an abrupt change in the bulk sound speed in the slab close to the Sunda Strait further suggesting a sharp change in the slab characteristics in this region. The possible causes of the observed rapid change in anisotropy, with reference to the other changes in geophysical characteristics are discussed further below. in geometry The Indian and Australian plates are subducting northwards beneath the Java-Sunda arc. As a result subduction moves from normal to the trench beneath Java to oblique subduction beneath Sumatra (most obviously characterised by the Great Sumatran transform fault). This change in the geometry of subduction across the region could potentially explain the change in anisotropy observed in the splitting results.

40- 80

 80 Ma and beneath Java 80-140 Ma. Other features of the slab are also different between these regions. Seismicity beneath Sumatra only extends to a depth of ∼300 km and as the slab bends beneath Java it extends to the base of the transition zone and this change in seismicity is located approximately at the same location as an abrupt change in the bulk sound speed (Gorbatov Page 21 of 47 A c c e p t e d M a n u s c r i p t and Kennett, 2003) and change in the location of the volcanic front (Syracuse

Page

  Institute, University of Tokyo. This research was funded by a Japan Society for the Promotion of Science (JSPS) postdoctoral fellowship (short term), JSPS/FF1/367 s c r i p t parameters in the presence of two anisotropic layers. Geophys. J. Int. 119, 949-963. Smith, G., Wiens, D., Fischer, K., Dorman, L., Webb, S., Hildebrand, J., 2001. A complex pattern of mantle flow in the Lau backarc. Science 292 (5517), 713-716. Syracuse, E., Abers, G., 2006. Global compilation of variations in slab depth beneath arc volcanoes and implications. Geochem. Geophys. Geosyst. 7. Teanby, N. A., Kendall, J.-M., Van der Baan, M., April 2004. Automation of shear-wave splitting measurements using cluster analysis. Bull. Seis. Soc. Am. 94 (2), 453-463. Wheller, G., Varne, R., Foden, J., Abbott, M., 1987. Geochemistry of Quaternary volcanism in the Sunda-Banda arc, Indonesia, and three-component genesis of island-arc basaltic magmas 32 (1-3), 137-160. Wolfe, C. J., Solomon, S. C., 1998. Shear-wave splitting and implications for mantle flow beneath the MELT region of the East Pacific Rise. Science 280, 1230-1232.

Figure 1 :

 1 Figure1: Map showing stations and events used in this study. White inverted triangles mark stations, blue circles show earthquakes (local events on main map, teleseismic events on top right map). Also shown are quaternary volcanos (red triangle), absolute plate motions[START_REF] Gripp | Current plate velocities relative to the hotspots incorporating the Nuvel-1 global plate motion model[END_REF], black arrows), slab contours at 100 km depth

Figure 2 :

 2 Figure 2: Example of the Teanby et al. (2004) method of shear-wave splitting. Top) SKS-wave splitting result at PSI and, bottom) local S-wave splitting result for PSI. a, g) Original traces (E, N, Z). b, h) Traces rotated into R and T directions before and after the anisotropy correction. R component is the initial shear-wave polarization before entering the anisotropic region. T component is perpendicular to the R component.Energy on the corrected transverse trace should be minimised in the analysis window. c, i) The top traces show the fast/slow shear waveforms for uncorrected (left) and corrected (middle(normalised)/right(real amplitudes)) seismograms. The bottom panels show the particle motion for uncorrected (left) and corrected (right) seismograms. A good result will show similar fast/slow shear waveforms and any elliptical particle motion will have been linearised. d, j) The results of the grid search over δt and φ. The method used minimises the second eigenvalue of the particle motion (i.e. the best result occurs where the particle motion is linear after removing the splitting). The optimum splitting parameters are represented by the cross, and the 1st surrounding contour denotes the 95% confidence region. e, k) Measurements of δt and φ obtained from 100 different analysis windows plotted against window number. A stable result will have a similar solution independent of analysis window. f, l) Cluster analysis of δt and φ obtained from 100 different analysis windows. A good result will have solutions clustered in one area.

Figure 3 :

 3 Figure 3: Local and SKS-wave splitting results beneath Sumatra (top), and Java (bottom). Yellow and blue bars show splitting result for SKS-wave arrivals and local S-waves respectively. Bars are oriented in the direction of the fast shear-wave polarisation, and the length of the bar is proportional to the amount of splitting. SKS-wave splitting results are plotted at the station location, local S-wave splitting results are plotted at the event location. Other symbols are as in Figure 1.

Figure 4 :

 4 Figure 4: Variation in local splitting delay time as a function of depth. Coloured crosses show local S-wave splitting delay times as a function of depth. SKS-wave splitting delay times are shown by coloured squares in the grey box at the bottom of the figure. Notethe lack of variation with depth for events <300 km deep, and the difference between delay time for events <300 km deep and events >300 km deep and SKS-events.

Figure 5 :

 5 Figure 5: a) Ray tracing through an isotropic slab model to find the path travelled from local events to the surface. Events are matched to the nearest point on a ray (local) or to their ak135 predicted piercing point (SKS/SKKS) b) The raypaths for stations UGM/SWH, showing the proportion of the ray spent in each part of the model. Splitting is predicted using the Christoffel equation( Červenỳ, 2001) for an appropriately rotated olivine tensor.

Figure 6 :

 6 Figure 6: Example shallow and slab anisotropy modelling result for PPI. For the shallow model anisotropy is confined to the uppermost 40 km of layer 1, representing the over-riding plate, and events from depths <300 km only are modelled. For the slab model anisotropy, the modelled anisotropy from the shallow anisotropy is fixed in layer 1, and a second layer of anisotropy exists in the slab. Events from depths >300 km and SKS/SKKS-phases are modelled. a,d) The calculated misfit surfaces for the grid search of anisotropy strength (i.e. degree of alignment -VFA) and a-axis orientation (γ) in a horizontal plane. Nodes in the grid are evaluated as the summed misfit across all event-station pairs, calculated for each δt and φ separately and also an equally weighted combined misfit. b,e) The comparison of modelled splitting (blue) and observed (black) for the optimum model, plotted at event location.c,f ) A density rotorgram of the misfit in a-axis direction for shallw and slab anisotropy respectively (γ; black is lowest misfit, white is highest), at the best fitting anisotropy strength VFA. This shows graphically the confidence in orientation of our best fitting model. Note the good fit to the data, and the well constrained direction.

Figure 7 :

 7 Figure 7: Variation in delay time as a function of (a) path length in the mantle wedge and over-riding plate (layer 1 in Figure 5) and b) path length above 410 km but beneath the

Figure 8 :

 8 Figure 8: Density rotorgrams (see Figure 6) for models of subduction zone anisotropy.

Figure 9 :

 9 Figure 9: A map showing the misfits (bars) for models with anisotropy in a thin layer above the slab (red), in the slab (green), and below the slab (blue), along with the best fitting anisotropy orientations (arrows). Note, the best fitting model places anisotropy in the slab with north-south orientations beneath Sumatra, and east-west orientations beneath Java.

  Figure 1:

  Figure 2:2

  Figure 4:

  Figure 8:

  Figure 9:

Table 1 :

 1 Best fitting models for shallow earthquakes (<300 km). Anisotropy is confined to a 40 km thick layer at the surface, and fits the data well with trench parallel anisotropy orientations.

Table 2 :

 2 Best fitting models for deep earthquakes (>300 km) and SKS phases. A model with anisotropy in the slab fits splitting results best. In all models the best result from the shallow models (Figure

Table 1 : 1 Table (

 11(

	s)
	Over-riding plate (layer 1) VFA χ 2 χ 2 /n -20 γ 0.20 46.92 6.70 -15 0.30 38.38 4.80 -55 0.40 1486.22 42.46 -20 0.20 26.98 5.40 -55 0.25 126.90 15.86 -20 0.50 268.79 67.20 -90 0.45 147.14 24.52 -35 0.25 1426.24 95.08 -80 0.30 1107.10 61.50 -75 0.25 95.14 47.57 2775.62 67.70 4769.81 44.16 events (n) γ VFA χ 2 χ 2 /n 11 0 1.00 43.79 3.98 2 30 0.65 248.01 124.00 11 15 0.55 531.06 48.29 6 85 1.00 269.93 44.99 4 0 1.00 57.74 14.435 2 5 0.45 315.19 157.60 36 1465.72 40.71 6 -90 1.00 295.74 49.29 7 -35 0.95 147.33 21.05 10 -15 0.85 296.16 29.616 23 739.23 32.14 59 2204.95 37.37 A c c e p t e d Number of Station events (n) BSI 7 GSI 8 PSI 35 PPI 5 KSI 8 KOTA 4 Sumatra misfits 67 LEM 6 UGM+SWH 15 BMI 18 KHK 2 Java misfits 41 Overall misfits 108 Supra-slab (layer 2) Number of 1994.19 29.76 γ 0 55 -30 10 5 10 -65 55 -75 M VFA 0.30 0.80 0.10 0.95 0.70 0.60 0.45 0.5 0.55 Slab (layer 3) Sub-slab (layer 4) χ 2 χ 2 /n γ VFA χ 2 44.88 4.08 0 0.1 48.01 3.92 1.96 55 0.3 5.20 568.75 51.70 -35 0.10 575.47 98.37 16.40 0 0.85 195.40 21.55 5.39 5 0.45 10.91 33.25 16.63 0 0.40 221.97 770.72 21.41 1056.96 81.89 13.65 75 0.5 125.68 114.37 16.34 -85 0.25 244.74 152.79 15.28 -75 0.25 146.57 349.05 15.18 516.99 1119.77 18.98 1573.95 a n u s c r i p t χ 2 /n 4.36 2.60 52.32 32.57 2.72 110.99 29.36 20.95 34.96 14.66 22.48 26.68
	Station	BSI	GSI	PSI	PPI	KSI	KOTA	Sumatra misfits	LEM	UGM+SWH	BMI	KHK	Java misfits	Overall misfits

Table 2 :

 2