Lev Truskinovsky

Anna Vainchtein

KINETICS OF MARTENSITIC PHASE TRANSITIONS: LATTICE MODEL

Keywords: Martensitic phase transitions, lattice models, nonlocal interactions, driving force, lattice waves, radiative damping AMS subject classifications. 37K60, 74N10, 74N20, 74H05

We obtain a closing kinetic relation for a mixed type hyperbolic-elliptic p-system originating in the theory of martensitic phase transitions by replacing continuum model with its natural discrete prototype. The procedure can be viewed as either regularization by discretization or as a physically motivated account of underlying discrete (atomic or mesoscopic) microstructure. Our fully inertial lattice model describes an isolated phase boundary and its novelty is in taking into account nonlocality in the form of general harmonic long-range interactions. Although the model is Hamiltonian at the microscale, it generates a nontrivial macroscopic jump relation between the velocity of the discontinuity and the conjugate configurational force. This relation characterizes the rate of (apparent) dissipation, respects entropy inequality but is supplementary to the usual Rankine-Hugoniot jump conditions. The dissipation at the macrolevel is due to the induced radiation of lattice waves carrying energy away from the propagating front. We show that sufficiently strong nonlocality has a significant effect on the kinetic relation in both near-sonic and small-velocity regions.

Introduction.

A characteristic feature of martensitic phase transitions in active materials is the energy dissipation leading to experimentally observed hysteresis. The dissipation is due to propagating domain boundaries which can be represented at the continuum level as surfaces of discontinuity. Classical elastodynamics admits nonzero dissipation on moving discontinuities but provides no information about its origin and kinetics. Although the arbitrariness of the rate of dissipation does not create problems in the case of classical shock waves, it is known to be the cause of nonuniqueness in the case of subsonic phase boundaries. The ambiguity at the macroscale reflects the failure of the continuum theory to describe phenomena inside the narrow transition fronts where dissipation actually takes place. The missing closing relation can be found by analyzing a regularized theory describing fine structure of the transition front. Usually the problem reduces to a study of a one-dimensional steady-state model (see [START_REF] Dafermos | Hyperbolic conservation laws in continuum physics[END_REF][START_REF] Lefloch | Hyperbolic systems of conservation laws[END_REF][START_REF] Serre | Systems of conservation laws[END_REF] for recent reviews).

To formulate the simplest problem of this type it is sufficient to consider longitudinal motions of a homogeneous elastic bar. The total energy of the bar is the sum of kinetic and potential contributions

E = ρ u2 2 + φ(u x ) dx, (1.1)
where u(x, t) is the displacement field, u ≡ ∂u/∂t is the velocity, u x ≡ ∂u/∂x is the strain, ρ is the mass density, and φ(u x ) is the elastic energy density. The function u(x, t) satisfies the nonlinear wave equation

ρü = (σ(u x )) x , (1.2)
where σ(u x ) = φ (u x ) is the stress-strain relation. Although in classical elastodynamics equation (1.2), also known as p-system, is hyperbolic, the hyperbolicity condition σ (u x ) > 0 is violated for martensitic materials with non-monotone stress-strain relation σ(u x ) [START_REF] Ericksen | Equilibrium of bars[END_REF]. This makes the initial-value problem associated with the mixed-type equation (1.2) ill-posed; in particular, it leads to the appearance of non-evolutionary or undercompressed discontinuities (kinks, phase boundaries, e.g. [START_REF] Ngan | Thermal trapping and kinetics of martensitic phase boundaries[END_REF][START_REF]Thermo-elastic aspects of dynamic nucleation[END_REF][START_REF]Shock Induced Transitions and Phase Structures in General Media[END_REF]). The parameters of these discontinuities satisfy both the classical Rankine-Hugoniot jump conditions

[[ u]] + V [[u x ]] = 0, ρV [[ u]] + [[σ(u x )]] = 0, (1.3)
where [[f ]] ≡ f +f -and V is the jump velocity, and the entropy inequality R = GV ≥ 0, where

G = [[φ]] -{σ}[[u x ]] (1.4)
is the configurational (driving) force and {σ} ≡ (σ + +σ -)/2, but fail to satisfy the Lax condition c + < V < c -, where c ± are the sound velocities in front and behind. One way to remedy the resulting instability is to supplement (1.3) by a relation specifying the dependence of the configurational force on the velocity of the phase boundary G = G(V ) (kinetic relation, [START_REF] Abeyaratne | A continuum model of a thermoelastic solid capable of undergoing phase transitions[END_REF][START_REF]Dynamics of nonequilibrium phase boundaries in a heat conducting elastic medium[END_REF]). Since the nonlinear wave equation ( 1.2) provides no information about the kinetic relation, the dependence G(V ) has often been either extracted from experiment [START_REF] Abeyaratne | On the kinetics of austenite→martensite phase transformation induced by impact in a Cu-Al-Ni shape-memory alloy[END_REF][START_REF] Escobar | On pressure-shear plate impact for studying the kinetics of stress-induced phase transformations[END_REF] or modeled phenomenologically [START_REF] Abeyaratne | A continuum model of a thermoelastic solid capable of undergoing phase transitions[END_REF][START_REF] Truskinovsky | Equilibrium interphase boundaries[END_REF][START_REF]Dynamics of nonequilibrium phase boundaries in a heat conducting elastic medium[END_REF]]. An alternative approach has been to derive the kinetic relation from an augmented model incorporating regularizing terms. A typical example is the viscosity-capillarity model, incorporating into (1.2) both dispersive and dissipative terms [START_REF] Slemrod | Admissibility criteria for propagating phase boundaries in a van der Waals fluid[END_REF][START_REF] Truskinovsky | Equilibrium interphase boundaries[END_REF]. The problem with this approach is that it introduces into the theory phenomenological dimensional parameters of unclear physical origin (see the review [START_REF] Fan | Dynamic flows with liquid/vapor phase transitions[END_REF]).

The aim of the present paper is to obtain the kinetic relation without any phenomenological assumptions at the macroscale by replacing the continuum model (1.2) with its natural discrete prototype. This procedure of going back from continuum to discrete can be viewed as either regularization by discretization or as a physically motivated account of underlying atomic or mesoscopic microstructure. While it is clear that the discrete model must be Hamiltonian to reproduce the conservative structure of the smooth solutions of (1.2), the energy dissipation on the discontinuities can be interpreted as the nonlinearity-induced radiation of lattice-scale waves which takes the energy away from the long-wave continuum level. This phenomenon is known in physics literature as radiative damping (e.g. [START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF]Lattice friction for crystalline defects: from dislocations to cracks[END_REF]).

Following this idea of regularization of (1.2) from "first principles", we consider fully inertial dynamics of a one-dimensional lattice combining bi-stability and longrange interactions. Following some previous work in fracture [START_REF] Slepyan | Dynamics of a crack in a lattice[END_REF] and plasticity [START_REF] Atkinson | Motion of a Frenkel-Kontorova dislocation in a onedimensional crystal[END_REF], we assume piecewise linear interactions allowing one to construct an explicit traveling wave solution of the discrete problem. There exists an extensive literature on shock waves and solitons in the local and nonlocal discrete systems with convex energy (e.g. [START_REF] Gaidideib | Effect of non-local interactions on soliton dynamics in anharmonic chains: Scale competition[END_REF][START_REF] Hou | Dispersive approximations in fluid dynamics[END_REF][START_REF] Peyrard | Discreteness effects on non-topological kink soliton dynamics in nonlinear lattices[END_REF][START_REF] Toda | Theory of nonlinear lattices[END_REF]) and on the semilinear analogs of the present system (e.g. Frenkel-Kontorova model [START_REF] Braun | Nonlinear dynamics of the Frenkel-Kontorova model[END_REF][START_REF]Oscillatory wave fronts in chains of coupled nonlinear oscillators[END_REF][START_REF] Peyrard | Simple theories of complex lattices[END_REF]). The discrete quasilinear problem for martensitic phase transitions in the lattices with nearest-neighbor (NN) interactions has been recently considered in [START_REF] Slepyan | Transition waves in bistable structures. II. Analytical solution: wave speed and energy dissipation[END_REF][START_REF]Explicit kinetic relation from "first principles[END_REF], following previous work on failure waves and shock waves in discrete structures [START_REF]Models and phenomena in Fracture Mechanics[END_REF][START_REF] Slepyan | Fracture wave in a chain structure[END_REF]. In the present paper we extend these results to the case of harmonic interactions of arbitrary range, including Coulomb or Kac-Baker interactions, among others. In addition to providing the formal solution to the general problem, we show that the local model with NN interactions (mean field approximation) is degenerate and consider in detail the simplest nonlocal model involving interaction between next-to-nearest neighbors (NNN). In particular, we study a transition from weak coupling case, when the domain of lattice trapping region is the largest, to an almost continuous, when the trapping region completely disappears. We also demonstrate that sufficiently strong NNN interactions generate a much broader class of admissible solutions than in the NN case with a possibility of radiation both in front and behind the moving discontinuity. The analytical solution obtained in this paper allows us to capture certain details that are difficult or impossible to detect in numerical simulations (e.g. [START_REF] Purohit | Dynamics of phase transitions in strings, beams and atomic chains[END_REF]), such as singular behavior of solutions near static-dynamic bifurcation, around resonances and in the sonic limit.

The paper is organized as follows. The discrete model with long-range interactions and the associated dynamical problem are introduced in Section 2. In Section 3 we formulate the dimensionless equations for the traveling waves, the boundary conditions and the admissibility conditions. An explicit solution for the traveling wave is obtained by Fourier transform in Section 4. In Section 5 we rederive known static solutions describing lattice-trapped states and link them to a nontrivial limit of the dynamic solutions. The energy transfer from long to short waves is studied in Section 6, where we obtain a closed-form kinetic relation. In Section 7 we illustrate the general theory on the case when the only long-range interactions are due to second nearest neighbors. The conclusions are presented in Section 8.

2. Discrete model. The simplest lattice structure can be modeled as a chain of particles connected with elastic springs. Suppose that the interactions are of longrange type and that every particle interacts with its q neighbors on each side. If u n (t) is the displacement of the nth particle, the total energy of the chain can be written as

E = ε ∞ n=-∞ ρ u2 n 2 + q p=1 pφ p u n+p -u n pε , (2.1)
where ε is the reference interparticle distance and φ p (w) is the energy density of the interaction between pth nearest neighbors. The dynamics of the chain is governed by an infinite system of ordinary differential equations:

ρü n = q p=1 φ p u n+p -u n pε -φ p u n -u n-p pε , (2.2)
which replace the nonlinear wave equation (1.2). The microscopically homogeneous configurations of the bar generate macroscopic stress-strain relation [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF] σ(w) = q p=1 pφ p (w).

(2.3)

For this correspondence to be valid, the functions φ p (w) must satisfy certain conditions preventing formation of microinhomogeneities which we specify for the case q = 2 in Section 7.

To obtain analytical results, we consider the simplest potentials compatible with phase transitions: bi-quadratic for local interactions (NN) and quadratic for nonlocal interactions (NNN, NNNN, etc.) Specifically we define

φ 1 (w) = 1 2 Ψ(1)w 2 , w≤ w c 1 2 Ψ(1)(w -a) 2 , w ≥ w c , (2.4)
and

φ p (w) = 1 2 pΨ(p)w 2 , p = 2, ..., q. (2.5)
The nonlinear springs representing NN interactions can be found in two different states depending on whether the strain w is below (phase I) or above (phase II) the critical value w c ; parameter a defines microscopic transformation strain. For simplicity we assume that the phases have equal elastic moduli Ψ(1) > 0. It is convenient to rewrite the governing equations (2.2) in terms of strain variables w n = (u nu n-1 )/ε and to use dimensionless variables:

t = t(Ψ(1)/ρ) 1/2 /ε w = w/a Ψ(p) = Ψ(p)/Ψ(1), p = 1, ..., q (2.6)
In terms of these variables with bars dropped, we can rewrite the governing equations (2.2) in the form ẅn -

|k-n|≤q Ψ(k -n)w k = 2θ(w n -w c ) -θ(w n+1 -w c ) -θ(w n-1 -w c ), (2.7)
where

Ψ(0) = -2 q p=1 Ψ(p), Ψ(-p) = Ψ(p) (2.8)
and θ(w) is a unit step function. The macroscopic stress-strain relation (2.3) takes the form

σ(w) = c 2 w -θ(w -w c ), (2.9)
where

c = q p=1 p 2 Ψ(p) 1/2 (2.10)
is the dimensionless macroscopic sonic speed. We choose the moduli Ψ(p) to ensure that the function σ(w) takes the form presented in Figure 2.1.

Traveling waves.

An isolated phase boundary moving with a constant velocity V can be obtained as a traveling wave solution of (2.7) with w n (t) = w(ξ), ξ = n -V t. We assume further that in the moving coordinate system all springs in the region ξ > 0 are in phase I (w n < w c ), and all springs with ξ < 0 are in phase II. The system (2.7) can then be replaced by a single nonlinear advance-delay differential equation:

V 2 w - |p|≤q Ψ(p)w(ξ + p) = 2θ(-ξ) -θ(-ξ -1) -θ(1 -ξ), (3.1)
The configurations at ξ = ±∞ must correspond to stable homogeneous equilibria plus superimposed short-wave oscillations with zero average:

w(ξ) → w ± as ξ → ±∞. (3.2)
The nonlinearity of the problem is in the switching condition

w(0) = w c . (3.3)
We assume that a solution is admissible if the springs in front of the moving interface remain in phase I and behind it in phase II. This implies that

w(ξ) < w c for ξ > 0, w(ξ) > w c for ξ < 0. (3.4)
Consequently, the mathematical problem reduces to solving (3.1) subject to (3.2), (3.3) and (3.4).

First observe that the equation (2.7) is linear in each phase (ξ > 0 and ξ > 0), which means that the solution can be represented as a superposition of linear waves w n = exp(i(knωt)). Since the elastic moduli are equal, the dispersion relation

ω 2 (k) = 4 q p=1 Ψ(p) sin 2 pk 2 (3.5)
is the same in both phases. In order for the linear modes to be compatible with the traveling wave ansatz, their phase velocity V p (k) = ω/k must be equal to V . This gives the restriction on the admissible wave lengths in the form

L(k, V ) = 0, (3.6)
where

L(k, V ) = 4 q p=1 Ψ(p) sin 2 pk 2 -V 2 k 2 . (3.7)
Among the modes selected by (3.6), the ones with complex wave numbers describe the core structure of the phase boundary, while the ones with real wave numbers correspond to radiation. Waves with k = 0 will be associated with the macroscopic configuration.

4. Exact solution. We solve Eq. (3.1) by writing w(ξ) = h(ξ) + w -and applying the complex Fourier transform ĥ

(k) = ∞ -∞ h(ξ)e i(k+iδ)ξ dξ, h(ξ) = 1 2π ∞+iδ -∞+iδ ĥ(k)e -ikξ dk,
where δ > 0 is a small parameter which guarantees convergence of the integrals. After inverting the Fourier transform and letting δ → 0, we obtain

w(ξ) = w -- 2 πi Γ sin 2 (k/2)e ikξ dk kL(k, V ) , (4.1)
where the contour Γ coincides with the real axis passing the singular point k = 0 from below. The singularities associated with nonzero real roots of L(k, V ) = 0 must comply with the radiation conditions. Specifically, the modes with group velocity V g = ∂ω/∂k larger than V can appear only in front, while the modes with V g < V can appear only behind the phase boundary [START_REF]Models and phenomena in Fracture Mechanics[END_REF]. Using the relation

V g = V + L k (k, V ) 2V k , (4.2)
where L k (k, V ) = ∂L/∂k and assuming V > 0, we obtain that V g ≷ V whenever kL k (k, V ) ≷ 0. Therefore to satisfy the radiation conditions, we need to dent the integration contour in (4.1) so that it passes below the singularities on the real axis

if kL k (k, V ) > 0 and above if kL k (k, V ) < 0.
To compute the integral (4.1) explicitly, we use the residue method closing the contour in the upper half-plane when ξ > 0 and in the lower half-plane when ξ < 0. The solutions look different in the generic case q > 1 and the degenerate case q = 1. For q > 1 the Jordan's Lemma can be applied directly and we obtain

w(ξ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ w -+ k∈M - 4 sin 2 (k/2)e ikξ kL k (k, V ) for ξ < 0 w -- 1 c 2 -V 2 - k∈M + 4 sin 2 (k/2)e ikξ kL k (k, V ) for ξ > 0, (4.3)
where

M ± = {k : L(k, V ) = 0, {Imk ≷ 0} {Imk = 0, kL k (k, V ) ≷ 0}}.
For q = 1 (short range interactions only) the contribution from a semi-arch at infinity does not vanish at ξ = ±0 and relations (4.3) must be supplemented by the following limiting conditions

w(ξ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ w -+ k∈M - 4 sin 2 (k/2)e ikξ kL k (k, V ) - 1 2 for ξ = 0- w -- 1 c 2 -V 2 - k∈M + 4 sin 2 (k/2)e ikξ kL k (k, V ) + 1 2 for ξ = 0+. (4.4)
In both cases, by applying the boundary conditions at infinity we obtain

w + = w -- 1 c 2 -V 2 . (4.5)
This condition coincides with the Rankine-Hugoniot relation

V 2 [[w]] = [[σ]],
for the macroscopic stress-strain relation (2.9). The continuity of w(ξ) at ξ = 0, implies that

1 c 2 -V 2 + k∈M 4 sin 2 (k/2) kL k (k, V ) = 1, q = 1 0, q > 1, (4.6)
where M = M + M -is the set of all nonzero roots of the dispersion relation (3.6). This condition is automatically satisfied for q > 1 since the sum of residues at all poles (including k = 0) equals zero; for q = 1 and ξ = 0 the integral over a contour at infinity contributes additional 1 in the right hand side of (4.6). The switching condition (3.3) together with (4.6) requires that

w ± = w c ∓ 1 2(c 2 -V 2 ) + k∈N pos 4 sin 2 (k/2) |kL k (k, V )| , (4.7)
where

N pos = N + pos N - pos , N ± pos = {k : L(k, V ) = 0, Imk = 0, k > 0, kL k (k, V ) ≷ 0}
are the sets of positive real roots with group velocities above and below V . By virtue of (4.5) the two conditions (4.7) are not independent and can be replaced by a single condition

1 2 (w -+ w + ) -w c = k∈N pos 4 sin 2 (k/2) |kL k (k, V )| . (4.8)
As we show in Section 6, equation (4.8) represents the desired kinetic relation. When V = 0, we can use the explicit formulae for w(ξ) to reconstruct the particle velocity profile v(ξ) = -V u (ξ). The relation between the velocity and the strain fields reads

v(ξ) -v(ξ -1) = -V w (ξ). (4.9)
where the right hand side is known. Solving (4.9) again by the Fourier transform, we obtain

v(ξ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ v + - V c 2 -V 2 -2V k∈M - sin(k/2)e ik(ξ+ 1 2 ) L k (k, V ) for ξ < - 1 2 v + + 2V k∈M + sin(k/2)e ik(ξ+ 1 2 ) L k (k, V ) for ξ > - 1 2 .
(4.10)

Notice that the average velocities at infinity satisfy the Rankine-Hugoniot condition (1.3) 1 , which in our case takes the form

v + -v -= V c 2 -V 2 .
To summarize, we have obtained a set of traveling wave solutions parametrized by the velocity V and the boundary value data w ± and v ± . The average particle velocity v + in front can always be set equal to zero due to the Gallilean invariance. If the strain in front of the discontinuity is also prescribed, the remaining three macroscopic parameters are fully constrained by the two classical Rankine-Hugoniot conditions plus a non-classical condition (4.8).

Static solutions.

A special treatment should be given to the case V = 0 when the continuous variable ξ = n -V t takes integer values, and the strain profile becomes discontinuous at every ξ = n. In this limit the differential equation reduces to a system of finite-difference equations, and we can replace the continuous Fourier transform by its discrete analog [START_REF] Celli | Motion of a screw dislocation in a crystal[END_REF][START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF]The relation between the solutions of mixed dynamical problems for a continuous elastic medium and a lattice[END_REF]. Observe that for a piecewise continuous function

ŵ(k) = ∞ -∞ w(x)e ikx dx = ∞ n=-∞ n+1 n w(x)e ikx dx.
Therefore, assuming that the strain profile w(x) converges to w 0 (n) as V → 0, we obtain

ŵ0 (k) = ∞ n=-∞ w 0 (n) e ik(n+1) -e ikn ik = e ik -1 ik ŵD 0 (k), (5.1)
where ŵD 0 (k) = ∞ n=-∞ w 0 (n)e ikn is the discrete Fourier transform of w 0 (n). Now we can use (4.1) to obtain

ŵ0 (k) = 2πδ(k)w -+ 4 sin 2 (k/2) ikω 2 (k) ,
where δ(k) is the Dirac delta function and ω 2 (k) is given by (3.5). Using (5.1) we can then find ŵD 0 (k) and applying inverse discrete Fourier transform, obtain a representation of the discrete solution:

w n = 1 2π π -π ŵD 0 (k)e -ikn dk = w -- 1 πi π -π
sin(k/2)e ik(n+1/2) dk ω 2 (k) .

To avoid the singularity at k = 0 we must pass it from below; all other roots of ω 2 (k) = 0 inside the strip -π ≤ Rek ≤ π have nonzero imaginary parts. Extending the integrand analytically by zero outside Rek ∈ [-π, π] and closing the contour of integration in the upper half-plane for n < 0 and lower half-plane for n > 0, we obtain by residue theorem

w n = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ w -- k∈F - sin(k/2)e ik(n+1/2) ω(k)ω (k) , n<0 w -- 1 c 2 + k∈F + sin(k/2)e ik(n+1/2) ω(k)ω (k) , n ≥ 0, (5.2)
where F ± = {k : ω 2 (k) = 0, Imk ≷ 0, -π ≤ Rek ≤ π}. Solutions satisfying the admissibility constraints

w n ≥ w c for n ≤ -1, w n ≤ w c for n ≥ 0 (5.3)
form a family of lattice-trapped equilibria (5.2) parametrized by w -, or, equivalently, by the total stress in the chain σ = c 2 w --1. The whole set of admissible stresses constitutes the trapping region. Notice that in this class of static solutions the phase boundary is pinned at the site n = -1. If the strain profile (5.2) is monotone, which occurs, for instance, when all long-range interactions are repulsive (Ψ(p) < 0 for p ≥ 2), the constraints (5.3) can be replaced by w 0 ≤ w c and w -1 ≥ w c . The trapping region can then be described explicitly

σ M -σ P ≤ σ ≤ σ M + σ P , (5.4)
where σ M = c 2 w c -1/2 is the Maxwell stress and

σ P = 1 2 -c 2 k∈F + sin(k/2)e ik/2 ω(k)ω (k) (5.5)
is the Peierls stress (see also [START_REF] Braun | Kinks in the Frenkel-Kontorova model with long-range interparticle interactions[END_REF][START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF]). The phase boundary remains trapped until the stress reaches one of the limiting values: σ = σ M -σ P , corresponding to w -1 = w c and interface moving to the left (V < 0) or σ = σ M + σ P , corresponding to w 0 = w c and interface moving to the right (V > 0). The two limiting solutions represent unstable equilibria from which the dynamic solution bifurcates. In Section 7 we compute the exact value of σ P for the case q = 2.

6. Kinetic relation. The waves generated in the core region carry the energy away from the front without changing the average values of parameters at infinity. At the continuum level these lattice waves are invisible and therefore the associated energy transfer is perceived as dissipation. To evaluate the rate of dissipation, we start with the microscopic energy balance

dE dt = A(t),
where E is the total energy of the chain and A(t) is the power supplied by the external loads. Since the solution of the discrete problem at infinity can be represented as a sum of the macroscopic contribution and the superimposed oscillations, we can split the averaged power accordingly. We obtain

A ≡ 1 τ τ 0 Adt = P -R, (6.1)
where P = σ + v +σ -v -is the macroscopic rate of work and R is the energy release due to radiated waves (interpreted at the macroscale as dissipation). While in the general case the expression for R may contain coupling terms, in the piecewise linear case the macroscopic and microscopic contributions decouple [START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF]Models and phenomena in Fracture Mechanics[END_REF]. Moreover, the dissipation rate R can be written as the sum of the contributions from the areas ahead and behind the front:

R(V ) = R + (V ) + R -(V ). (6.2)
Furthermore, due to the asymptotic orthogonality of the linear modes, the terms in the right hand side of (6.2) can be expressed in terms of contributions due to individual modes. Since the energy flux associated with the linear mode k is the product of the average energy density G k and the relative velocity |V g -V | of the energy transport with respect to the moving front, we can write

R + (V ) = k∈N + G k + (V g -V ), R -(V ) = k∈N - G k -(V -V g ). (6.3)
where

G k = G -G 0 k ≡ 1 τ τ (k) 0 (G -G 0 )dt with τ (k) = 2π/ω(k).
Here

G(ξ) = 1 2 (v(ξ)) 2 + c 2 2 (w(ξ)) 2 -θ(-ξ)(w(ξ) -w c ) - 1 2 q-1 p=1 B(p){(w(ξ + p) -w(ξ)) 2 + (w(ξ) -w(ξ -p)) 2 } (6.4)
with B(p) = q-|p| l=1 lΨ(l + |p|), and

G 0 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 2 v 2 -+ c 2 2 w 2 --w -+ w c , ξ < 0 1 2 v 2 + + c 2 2 w 2 + , ξ>0
is the corresponding energy density of the homogeneous states. To complete the computation we need to specify the asymptotic representation of the velocity and strain fields at ±∞. We can write

v(ξ) ≈ v 0 (ξ) + v k (ξ), w(ξ) ≈ w 0 (ξ) + w k (ξ).
where

v 0 (ξ) = v -, ξ < 0 v + , ξ > 0 , w 0 (ξ) = w -, ξ < 0 w + ξ > 0
are the homogeneous components and

v k (ξ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ - 4V sin(k/2) cos(k(ξ -1/2)) L k (k, V ) , ξ < 0, k ∈ N - pos 4V sin(k/2) cos(k(ξ -1/2)) L k (k, V ) , ξ >0, k ∈ N + pos , w k (ξ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 8 sin 2 (k/2) cos kξ kL k (k, V ) , ξ <0, k ∈ N - pos - 8 sin 2 (k/2) cos kξ kL k (k, V ) , ξ > 0, k ∈ N - pos (6.5)
are the oscillatory components. The average energy density carried by the wave with the wave number k ∈ N ± pos can now be written as

G k ± = 1 2τ (k) τ (k) 0 v 2 k (ξ) + c 2 (w k (ξ)) 2 - q-1 p=1 B(p){(w k (ξ + p) -w k (ξ)) 2 + (w k (ξ) -w k (ξ -p)) 2 )} dt = 8V 2 sin 2 (k/2) (L k (k, V )) 2 .
which gives for the total energy flux

R(V ) = k∈N + pos 4V sin 2 (k/2) kL k (k, V ) - k∈N - pos 4V sin 2 (k/2) kL k (k, V ) = k∈N pos 4V sin 2 (k/2) |kL k (k, V )| .
Recalling the definition R(V ) = G(V )V we obtain the microscopic expression for the configurational force:

G(V ) = k∈N pos 4 sin 2 (k/2) |kL k (k, V )| . (6.6)
The function G(V ) is now known since both L(k, V ) and N pos depend explicitly on V . Comparing (6.6) with the macroscopic definition of the configurational force (1.4) we obtain which can be interpreted geometrically as the area difference between two shaded triangles in Figure 2.1. Notice that by combining (6.6) and (6.7) we obtain exactly (4.8), which shows that (4.8) is indeed a kinetic relation and that micro and macro assessments of dissipation are compatible.

G = 1 2 (w -+ w + ) -w c , (6.7)

Example.

To illustrate the general solution obtained above, in what follows we consider a special case q = 2 when both NN and NNN interaction are taken into account (see Figure 7.1). The model is then characterized by two elastic constants Ψ(1) and Ψ(2). It can be shown [START_REF]Quasi-continuum modeling of short-wave instabilities in crystal lattices[END_REF] that the homogeneous equilibria in this model are stable if and only if

Ψ(1) > 0, Ψ(1) + 4Ψ(2) > 0. (7.1)
We assume that these conditions are satisfied and rescale the equations by introducing a single dimensionless parameter

β = 4Ψ(2)/Ψ(1),
measuring the relative strength of NNN and NN interactions. The stability constraints give -1 < β ≤ ∞; one can further restrict this interval to

-1 < β ≤ 0 (7.2)
if we recall that the inequality Ψ(2) < 0 is suggested by the linearization of the potentials of the Lennard-Jones type [START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF][START_REF]The origin of nucleation peak in transformational plasticity[END_REF]. The total energy of the system now reduces to

E = ∞ n=-∞ v 2 n 2 + 1 + β 2 w 2 n -θ(w n -w c )(w n -w c ) - β 4 (w n+1 -w n ) 2 . (7.3)
According to (7.3), parameter β/(1 + β) characterizes the effect of discreteness: if β ∼ 0 we have the case of weak coupling (essential discreteness) while if β ∼ -1 we have the case of strong coupling implying quasicontinuum behavior. This is compatible with the fact that at β = 0 the Peierls stress characterizing the width of the latticetrapping domain takes the largest value (spinodal limit), while at β = -1 the Peierls stress is zero and the trapping region disappears. The energy 7.3 produces the following equation for the traveling waves

V 2 w - β 4 w(x + 2) -2w(x) + w(x -2) -w(x + 1) + 2w(x) -w(x -1) = 2θ(-x) -θ(-x -1) -θ(1 -x). (7.4)
The formal solution of this equation has been obtained in Section 4. Below we make it explicit and provide illustrations for the physically justified range of parameters.

Dispersion relation.

To compute the strain and velocity profiles at a given V we need to find the nonzero roots of the dispersion relation

L(k, V ) = 4 sin 2 (k/2) + β sin 2 k -V 2 k 2 = 0. (7.5)
We denote the roots by k = k 1 + ik 2 and divide them into three categories: real, responsible for radiation; purely imaginary, providing the monotone structure of the interphase region and complex, describing oscillatory contributions to the core. Since L(k, V ) is an even function of k, the real roots appear in pairs k = ±k 1 . Assuming positive V , we obtain

V (k 1 ) = 4 sin 2 (k 1 /2) + β sin 2 k 1 |k 1 | .
This function is plotted in Figure 7.2a. An infinite number of local maxima on this graph V = V i , correspond to resonance velocities: at these points L k (k, V ) = 0 and the sums in (4.3), (4.7) and (4.10) diverge. Between the resonance velocities, equation (7.5) possesses a finite number of positive real roots corresponding to propagating waves. To determine whether these waves propagate ahead or behind the front, we need to check whether kL k (k, V ) = 2k 3 V (k)V (k) is positive or negative. At V > 0 the radiation conditions say that the waves with kV (k) > 0 propagate in front of the phase boundary, while the waves with kV (k) < 0 propagate behind. The value of β affects the function V (k) noticeably only at long waves (small k). Observe that V (0) is equal to the macroscopic sound speed c = (1 + β) 1/2 , that V (0) = 0 and that

V (0) = - 1 + 4β 12 √ 1 + β .
At -1/4 < β ≤ 0, the function V (k) has a maximum at k = 0 while at -1 < β < -1/4 it has a local minimum implying that sufficiently strong coupling (β < -1/4) creates the possibility of the microscopic waves waves moving faster than the macroscopic sound speed. The range of supersonic speeds increases as β → -1, and in the limiting case β = -1 all propagating waves are macroscopically supersonic. It is interesting that the critical value β = -1/4 also emerges in the strain-gradient approximation of the energy (7.

3), where it corresponds to the change of sign of the coefficient in front of the higher-order term [START_REF] Mindlin | Second gradient of strain and surface tension in linear elasticity[END_REF][START_REF] Triantafyllidis | The influence of scale size on the stability of periodic solids and the role of associated higher order gradient continuum models[END_REF]. In this approximation the dispersion relation V (k) is replaced by a parabola: for β > -1/4 (weak nonlocality) the parabola is directed downward and the strain-gradient coefficient is negative while for β < -1/4 (strong nonlocality) the parabola is upward and the strain-gradient contribution to the energy is positive definite. The latter implies that subsonic waves are absent, and in order to yield a nontrivial kinetic relation the corresponding quasicontinuum model must be augmented by even higher-order terms [START_REF]Kinetics of lattice phase transitions II: Quasicontinuum approximation[END_REF].

The purely imaginary roots of (7.5) also appear in symmetric pairs and correspond to non-oscillatory modes exponentially decreasing away from the front. By solving L(ik 2 , V ) = 0 for V (k 2 ) we obtain

V (k 2 ) = 4 sinh 2 (k 2 /2) + β sinh 2 k 2 |k 2 |
This function is shown in Figure 7.2b. A straightforward computation shows that

V (0) = c, V (k 2 )| k 2 =0 = 0 and V (k 2 )| k 2 =0 = (1 + 4β)/(12 √ 1 + β). For -1 < β < -1/4 the maximum of the curve V (k 2
) is reached at k = 0, which means that in the case of strong coupling only macroscopically subsonic phase boundaries have monotone contribution to the core structure. Both the value V (0) and the range of available wave numbers decrease as β → -1, so that in the limit purely imaginary roots disappear. In the case β = 0 the function V (k 2 ) is convex everywhere, implying that no imaginary roots contribute to the subsonic solution. This is compatible with the fact that in the degenerate NN limit the static interface (V = 0) becomes atomically sharp.

The complex roots contribute to the oscillatory structure of the core region. For the given V the real (k 1 ) and imaginary (k 2 ) parts of the dispersion spectrum satisfy the system of two equations: ReL(k, V ) = 0 and ImL(k, V ) = 0. The set of complex roots contains infinitely many branches that come in symmetric quadruples. The first quadrant of the complex plane is shown in Figure 7.3. Observe that the complex roots can be divided into two sets: Q and P . The set Q (thick gray lines), has a purely dynamic nature, exists for all β and consists of branches bifurcating from the local maxima of the real branch and (for β < -1/4) from an imaginary branch. These roots contribute to the boundary layers around the front only at nonzero V . The second set of complex roots, P , shown in Figure 7.3 by thick black lines contains purely imaginary branches and bifurcates of complex roots intersecting the plane V = 0. These branches contribute to the static solution: at V = 0 they can be given explicitly given by

k = 2πn ± iλ, λ = 2arccosh 1 |β| , (7.6)
where n is an integer [START_REF]The origin of nucleation peak in transformational plasticity[END_REF]. As β tends to zero, the imaginary parts of P -roots approach ±∞; the eventual disappearance of these roots in the limit β → 0 is responsible for the sharpening of the front in the NN approximation.

Strain and velocity profiles.

Typical profiles of strain w(x) and velocity v(x) computed for the NNN model from (4.3), (4.10) are shown in Figure 7.4 where β = -0.2 and the first two resonance velocities are V 1 = 0.2164 and V 2 = 0.1282. At V = 0.5 > V 1 we obtain only one radiative mode propagating behind the phase boundary; at V 2 < V = 0.16 < V 1 the solution exhibits two additional radiative modes, one propagating behind and one in front of the phase boundary.

A closer inspection of the solutions at V < 0.266 reveals a violation of the constraints (3.4). For example, in the strain profile corresponding to V = 0.16 in Figure 7.4 the maximum strain w > w c is achieved at some point ξ > 0. Our numerical computations suggest that the entire velocity interval (0, 0.266) around the resonances may have to be excluded (at this particular β). Similar "velocity gaps" were also detected in [START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF]Lattice friction for crystalline defects: from dislocations to cracks[END_REF][START_REF] Marder | Origin of crack tip instabilities[END_REF] for the semilinear Frenkel-Kontorova problem.

At larger β velocity gaps become narrower and steady interface propagation becomes possible in certain subcritical velocity intervals. For example, at β = -0.75 traveling wave solutions exist in the intervals: [0.24, 0.5] (between the first resonance V 1 = 0.2150 and the sonic speed c = 0.5); [0.142, 0.19] (between first and second resonances, V 2 = 0.1279); [0.1, 0.11] (between V 2 and V 3 = 0.0912); [0.078, 0.08] (between V 3 and V 4 = 0.0708) and possibly in some shorter intervals at smaller V . Two such solutions are shown in Figure 7.5. The first one corresponds to V = 0.16 which is between the first and second resonances; unlike its counterpart at β = -0.2, this solution is now admissible. The second admissible profile corresponds to the value of velocity V = 0.105 located between second and third resonances. In this case there are five radiative modes, two in front and three behind the phase boundary.

The growth of small-velocity intervals where steady interface propagation exists is due to the presence of P -roots at nonzero β. As β becomes larger in absolute value, these roots move closer to the real axis, widening the transition layer and suppressing oscillations due to other complex roots. The effect of varying β on existence of smallvelocity solutions may thus be compared to the effect of viscosity in the Frenkel-Kontorova model [START_REF] Caprio | Depinning transitions in discrete reaction-diffusion equations[END_REF][START_REF] Kresse | Lattice models of propagating defects[END_REF]: both stronger nonlocality and larger viscosity reduce the velocity gaps. 7.3. Kinetic relation. Using (6.6) and the known dispersion spectrum, we can now explicitly evaluate the kinetic relation G(V ). A representative example is shown in Figure 7.6. At resonance velocities the configurational force required to move the interface tends to infinity which explains the appearance of the peaks. As we discussed above, at sufficiently small β the entire region around the small-velocity resonances has to be excluded since the corresponding solutions violate the constraints (3.4); see Figure 7.7a. As β increases, some of the small-velocity solutions between the resonances become admissible, as shown in Figure 7.7b-d.

Zero-velocity limit. To obtain solution at V = 0 we can use (5.2) with F ± = {±iλ}, where λ is defined in (7.6). After some algebraic manipulations, the family of latticetrapped equilibria (5.2) can be represented in the form

w n = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ σ + 1 1 + β - e λ(n+1/2) 2(1 + β) cosh(λ/2) , n < 0 σ 1 + β + e -λ(n+1/2) 2(1 + β) cosh(λ/2) , n ≥ 0, (7.7)
where σ lies in the region (5.4). The solutions (7.7) coincide with equilibria obtained by a different method in [START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF] where it was shown that they correspond to local minima of the energy and are therefore metastable. The expression for the Peierls stress (5.5) indicating the boundary of the metastability region can now be written explicitly as

σ P = 1 2 1 + β.
Observe that at β = 0 the Peierls stress coincides with the spinodal stress σ S = 1/2. As β grows, the trapping region becomes narrower, and eventually disappears at β = -1 (Peierls stress is zero). The upper boundary of the trapping region (5.4) corresponds to the case when w 0 = w c , which coincides with condition (3.3) satisfied by the dynamic solutions at V > 0. This saddle-point configuration can be given explicitly

w n = lim V →0 w(n -V t) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ w c + e λ/2 -e λ(n+1/2) 2(1 + β) cosh(λ/2)
, n < 0

w c + e -λ(n+1/2) -e -λ/2 2(1 + β) cosh(λ/2) , n ≥ 0. (7.8)
Using this solution, we can also obtain the limiting value of the configurational force corresponding to the Peierls depinning limit (static-dynamic bifurcation)

G(0) = 1 2 (w -+ w + ) -w c = 1 2 √ 1 + β . (7.9)
Observe that although the Peierls stress tends to zero when β → -1, the corresponding value of the configurational force G(0) tends to infinity, which is an artifact of our oversimplified model and is due to the divergence of the macroscopic transformational strain in this limit.

Sonic limit. The qualitative behavior of the function G(V ) in the limit V → c depends on β. If -1/4 < β ≤ 0, and V c, the wave spectrum contains a single wave number k which approaches zero as V → c. Expanding the expression for the configurational force (6.6) at small k we obtain Figures 7.6 and 7.7a). The picture is qualitatively different when -1 < β < -1/4. In this case as V approaches c from below, the limit of the corresponding real wave number k is nonzero k s , and therefore configurational force G(V ) remains finite (see Figure 7.7b-d).

G ∼ 6 (1 + 4β)k 2 , which implies that G(V ) → ∞ as V → c (see
Going back to the general case of arbitrary V < c, we observe that there is an infinite number of critical values of β at which the sonic speed c coincides with one of the resonance velocity. For example, at β = -0.9539 we have c = V 1 = 0.2147 implying that for β ≤ -0.9539 the subsonic region lies below the first resonance (see Figure 7.6f). Although the total domain of existence of the traveling wave solutions shrinks as β → -1 the domain of admissible traveling waves between the resonances expands. Our numerical calculations for the case q = 3 (NNNN model) show qualitatively similar behavior of the kinetic curves, suggesting that the present model may be already capturing the main effects of nonlocality.

Conclusions.

We used a simple discrete model to derive an explicit kinetic relation in one-dimensional theory of martensitic phase transitions. The macroscopic dissipation was interpreted as the energy of the lattice waves emitted by a moving interface, and we showed that discrete and continuum assessments of the dissipation are compatible. The present study complements previous analyses of the local model with NN interactions by including an arbitrary number of harmonic long-range interactions.

We showed that the local model is degenerate and analyzed in detail the nonlocal model accounting for next-to-nearest neighbor (NNN) interaction. By varying the strength of the NNN interactions we could pass from the essentially discrete case to the almost continuum case. We showed that nonlocality affects the size of lattice trapping: as NNN bonds become stronger, the trapping region reduces in size in terms of stresses. However, it increases in terms of driving forces, which emphasizes an important difference between the real and configurational forces. In addition to enlarging the domain of existence of steady state regimes with high-frequency radiation in both directions, sufficiently strong coupling significantly alters the mobility curves near the sonic speeds. Contrary to the simplest local theory, strongly nonlocal model produces multivalued kinetic relations with several admissible branches and rich variety of configurations of emitted lattice waves. Although we have not studied stability of the constructed traveling waves, the fact that the dissipation potential associated with these dynamic regimes is locally convex suggests that they may be infinitesimally stable.

The present work was motivated by recent studies of crack and dislocation dynamics (e.g. [START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF]Lattice friction for crystalline defects: from dislocations to cracks[END_REF]), where a two-well version of the discrete Frenkel-Kontorova (FK) model was used to compute the macroscopic rate of dissipation. The formal difference between the semilinear FK model and the quasilinear model presented in this paper is that now nonlinearity concerns discrete derivatives. While the two problems are equivalent in statics and in overdamped dynamic limit [START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF]40], the fully inertial models are different. An additional level of complexity in the quasilinear setting is associated with the presence of the limiting characteristic velocity, nonzero macroscopic particle velocity and the necessity to satisfy the discrete analogs of the Rankine-Hugoniot jump conditions.

2 A 1 Fig. 2 . 1 .

 2121 Fig. 2.1. The bi-linear macroscopic stress-strain law and the Rayleigh line connecting the states at infinity for a martensitic phase boundary. The difference between the shaded areas A 2 -A 1 represents the configurational force.

Fig. 7 . 1 .

 71 Fig. 7.1. One-dimensional discrete microstructure with nearest and next-to-nearest-neighbor interactions (q = 2).

Fig. 7 . 2 .

 72 Fig. 7.2. Real (a) and imaginary (b) roots of the dispersion relation L(k, V ) = 0 at different β.

1 Fig. 7 . 3 .

 173 Fig. 7.3. The structure of nonzero roots of L(k, V ) = 0 in the cases of strong (β = -2/3) and weak (β = -1/16) nonlocality. Thin lines -real roots, set N ; thick lines -P -roots; gray lines -Q-roots.

Fig. 7 . 4 .

 74 Fig. 7.4. Strain and velocity profiles at β = -0.2, w c = 1, v + = 0.

Fig. 7 . 5 .Fig. 7 . 6 .

 7576 Fig. 7.5. Strain profiles at β = -0.75, w c = 1 and V < V 1 .

Fig. 7 . 7 .

 77 Fig. 7.7. Kinetic relations G(V ) for different β.

Acknowledgements. This work was supported by the NSF grants DMS-0102841 (L.T.) and DMS-0137634 (A.V.).


Marc Bonnet -Bojan

B Guzina
email: guzina@wave.ce.umn.edu

Topological derivative applied to cavity identification from elastodynamic surface measurements

Keywords: cavity identification, elastodynamics, inverse scattering, adjoint field method, shape sensitivity, topological derivative. MOTS-CLÉS : identification de cavité, élastodynamique, diffraction inverse, méthode de l'état adjoint, dérivée géométrique, dérivée topologique

published or not. The documents may come

Topological derivative applied to cavity identification from elastodynamic surface measurements

Introduction

Three-dimensional imaging of cavities embedded in a solid using elastic waves is a topic of intrinsic interest in a number of applications ranging from nondestructive material testing to oil prospecting and underground object detection. In situations when detailed mapping of objects (defense facilities, buried waste) is required and only a few measurements can be made, the use of surface discretization-based boundary integral equation (BIE) techniques provides the most direct link between the surface measurements and the buried geometrical objects. Such an approach is well established for acoustic problems [START_REF] Kress R | Inverse acoustic and electromagnetic scattering theory[END_REF]. More recently, a BIE-based analytical and computational framework for the identification of cavities in a semi-infinite solid from surface elastodynamic measurements has been developed and reported in [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF][START_REF] Nintcheu | A Computational Basis for Elastodynamic Cavity Identification in a Semi-Infinite Solid[END_REF]. Computing time considerations are such that global search methods are currently impractical for three-dimensional elastodynamic inverse scattering problems, and hence led to using a conventional gradient-based 'blind' optimization scheme. As often with such schemes, the highly non-convex character of the cost function makes the identification result dependent on the initial guess.

This has prompted the authors to investigate the use of topological derivative as a tool for preliminary probing. The formulation presented in this article is applicable to bounded or unbounded elastic media. A set of numerical results is included to illustrate the utility of topological derivative for outlining the cavity location, shape and size prior to doing an actual inversion of measurements. Despite the fact that the results presented here were obtained from a BIE solution, the proposed methodology is readily applicable to other computational platforms such as the finite element method.

Preliminaries

Let Ω true denote an elastic body bounded by the external surface S, divided into complementary subsets S N and S D supporting prescribed tractions and displacements, respectively. An unknown cavity (or a set thereof) B true bounded by the closed surface(s) Γ true is embedded in Ω true , so that Ω true = Ω\ Btrue where Ω is the reference, i.e. cavity-free, counterpart of Ω true . On applying a steady-state traction p on S N with angular frequency ω, an elastodynamic state u true arises which solves the problem L(ω)u true = 0 (in Ω true ), t true = p (on S N ), t true = 0 (on Γ true ), u true = 0 (on S D ).

[1]

Here ρ is the mass density, C is the fourth-order isotropic elasticity tensor characterized by the shear modulus μ and Poisson's ratio ν, L(ω) is the Navier linear partial differential operator, i.e. L(ω)w = div(C : ∇w) + ρω 2 w, and t true ≡ σ true •n = (C : ∇u true )•n denotes the traction vector associated with the displacement u true through Hooke's law. For simplicity, it is assumed that ω is not an eigenfrequency of any of the boundary-value problems appearing in the ensuing developments.

For the inverse problem of interest, where the location, topology and geometry of B true (or equivalently Γ true ) is being sought, the trace of u true on S (denoted hereafter as u obs ) is assumed to be available over the measurement region S obs ⊂ S N . Let u c denote the solution to the forward problem for a given excitation p and a trial cavity B c bounded by Γ: u c is then defined over Ω c = Ω\ Bc , and governed by the equations

L(ω)u c = 0 (in Ω c ), t c = p (on S N ), t c = 0 (on Γ). u c = 0 (on S D ). [2]
where t c is the traction vector associated with u c . To solve the inverse problem, a misfit cost function is set up in order to minimize the difference between u obs and u c . Generic cost functions having the format

J (Ω c ) = S obs ϕ(u c (ξ), ξ) dΓ ξ [3]
are considered. The output least-squares cost function, commonly used for such purpose, corresponds to the particular case where 2ϕ(w, ξ)

= (w(ξ)-u obs (ξ)) • W (ξ) • (w(ξ)-u obs (ξ))
, where W (ξ) is a 3 × 3 matrix-valued weighting function, assumed to be symmetric and positive definite (the simplest choice being W (ξ) = I 2 ), while overbar denotes complex conjugation.

Topological derivative

To aid the gradient-based minimization of J (Ω), a tool often used for identifying B true on the basis of u obs , the development of topological derivative for the cost functionals of form [3], which would facilitate a rational selection of the necessary initial "guess" in terms of the location, topology and geometry of B true , is investigated. To this end, let B ε (x o ) = x o + εB, where B ⊂ R 3 is a fixed bounded open set with boundary S and volume |B| containing the origin, define the region of space occupied by a cavity of (small) size ε > 0 containing a fixed sampling point x o . Following [GAR01, SOK99], one is in particular interested in the asymptotic behavior of J (Ω ε ) for infinitesimal ε > 0, where Ω ε = Ω\B ε (x o ), and

B ε (x o ) is the closure of B ε (x o ).
The topological derivative T (x o ) of the cost functional J (Ω) at x o for a cavity-free body is hence defined through an expansion of the form:

J (Ω ε ) = J (Ω)+f (ε) |B| T (x o )+o(f (ε)) (ε Diam(Ω), B ε (x o ) ⊂ Ω) [4]
where the function f defines the leading asymptotic behavior of J (Ω ε ) and is such that f (ε) → 0 as ε → 0. The definition [4] is not restricted to spherical infinitesimal cavities (for which B is the unit ball, S the unit sphere and |B| = 4π/3). In general, the value T (x o ) is expected to depend on the shape of B.

With reference to [4], the evaluation of J (Ω ε ) requires the knowledge of the elastodynamic solution u ε to the forward problem [2] with B c replaced by

B ε ≡ B ε (x o ).
It is thus convenient to decompose the total displacement field u ε as u ε = u + ũε , where ũε denotes the scattered field and u is the free field defined as the response of the void-free (reference) solid Ω due to given excitation (i.e. boundary traction) p, so that

L(ω)u = 0 (in Ω), t = p (on S N ), u = 0 (on S D ), [ 5
]
where t is the traction vector associated with u, and

L(ω)ũ ε = 0 (in Ω ε ), tε = 0 (on S N ), ũε = 0 (on S D ), tε = -(C : ∇u)•n (on Γ ε ) [6]
where σ = C : ∇u is the stress tensor associated with the free field [5], Γ ε is the boundary of B ε , and n is the normal on S ∪ Γ ε outward to Ω ε . For infinitesimal ε the scattered field is expected to vanish, i.e. lim ε→0 |ũ ε (x)| = 0 (x ∈ Ω c ), whereas the free-field, by its definition [5], does not depend on ε. One may expand J (Ω ε ) with respect to ũε as

J (Ω ε ) = S obs ϕ(u ε (ξ), ξ) dΓ ξ = S obs ϕ(u(ξ), ξ) + Re ∂ϕ ∂u u(ξ), ξ • ũε (ξ) + o(|ũ ε (ξ)|) dΓ ξ = J (Ω) + S obs Re ∂ϕ ∂u u(ξ), ξ • ũε (ξ) dΓ ξ + o( ũε ) [7]
where

∂ϕ ∂w ≡ ∂ϕ ∂w R -i ∂ϕ ∂w I w R = Re(w) , w I = Im(w) [8]
By means of [4] and [7], the topological derivative of J (Ω) can be recast as:

T (x o ) = lim ε→0 1 f (ε) |B| S obs Re ∂ϕ ∂u u(ξ), ξ • ũε (ξ) dΓ ξ . [9]
One is then left with the task of evaluating the leading contribution of the integral in [9] for ε 1.

Adjoint field approach

Define the adjoint field u(ξ) as the solution to the elastodynamic problem

L(ω) u = 0 (in Ω), t = ∂ϕ ∂u (u, •) (on S obs ), t = 0 (on S N \S obs ), u = 0 (on S D ). [10]
where u is the free-field defined by [5], the prescribed traction is defined in terms of the cost function density ϕ, and the convention [8] is employed. Application of the elastodynamic reciprocity theorem to u and ũε over Ω ε yields

Γ ε ( u• tε -ũε • t) dΓ ξ = § obs ũε • ∂ϕ ∂u (u, •) dΓ ξ , [11]
which makes use of the boundary conditions in [6] and [10]. In practical terms, the assumed distribution of p over S obs can be interpreted as being proportional to a measure of the misfit between the measured displacement field u obs and the (reference) free-field u. In particular, in the least-squares case, ∂ϕ/∂u = W T (uu obs ) so that p depends linearly on the observed scattered field uu obs . On the basis of [6], [7] and [11], one finds that

J (Ω ε ) = J (Ω) + Re Γ ε u • tε dΓ - Γ ε ũε • t dΓ + o( ũε ) = J (Ω) -Re Γ ε u •t dΓ + Γ ε ũε • t dΓ + o( ũε ) as ε → 0. [12]
For a given primary excitation p, the first integral over Γ ε in [12] is known since the adjoint field u is determined solely in terms of u, u obs and the cost function density ϕ (see [10]). In fact, by virtue of [5], one has

Γ ε u•t dΓ = - B ε ∇•( u•C : ∇u) dΩ = B ε ρω 2 u•u -∇ u : C : ∇u dΩ = ε 3 |B| ρω 2 u•u -∇ u : C : ∇u (x o ) + o(ε 3 ) as ε → 0 [13]
where the minus sign in front of the first integral over B ε appears because t is defined in terms of the unit normal n pointing to the interior of B ε .

To complete the expansion of

J (Ω ε ) in [12]
, one is left with evaluating the leading asymptotic behavior of the last integral for vanishing ε. To accomplish this task, the asymptotic behavior of ũε on Γ ε must be investigated. This study [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] is based on an asymptotic analysis of the governing BIE for ũε and yields the following result:

ũε (ξ) = εσ k (x o )U k ( ξ) + o(ε) (ξ ∈ Γ ε ) [14] with ∇ξ •(C : ∇ξU k ) = 0 ( ξ ∈ R 3 \ B), (C : ∇ξU k )•n = -(n k e + n e k )/2 ( ξ ∈ S ),
where ξ = ξ/ε is a normalized position vector. The six canonical displacements U k ( ξ) are seen to solve normalized exterior elastostatic problems that do not depend on x o and ε. Then, by employing [14], one obtains

Γ ε ũε • t dΓ = ε 3 σ ij (x o ) σ k (x o ) S U k i ( ξ) n j ( ξ) dϑξ +o(ε 3 ) as ε → 0. [15]
where dϑξ denotes the area differential element on the unit sphere. On substituting [15] and [13] into [12], the leading asymptotic behavior of J (Ω ε ) for vanishing ε is found to take the form

J (Ω ε ) = J (Ω) + ε 3 |B| Re σ : A : σ -ρω 2 u•u (x o ) . [16]
Note in particular, with reference to [4], that the analysis shows that f (ε) = ε 3 , which is also known to hold for traction-free cavities in the 3-D elastostatic case [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF][START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF]. The constant fourth-order tensor A is defined by

A ijk = 1 2μ I ijk - ν 1 + ν δ ij δ k - 1 |B| S U k i ( ξ)n j ( ξ) dϑξ [17]
Finally, by virtue of [4] and [16], the sought adjoint-field formula for the topological derivative is

T (x o ) = Re σ : A : σ -ρω 2 u•u (x o ) [18]
For an arbitrarily shaped infinitesimal cavity, the six canonical problems in [14] should in general be solved numerically. This is a modest computational task, since in fact one only needs to solve six elementary static problems which do not depend on x o . For the particular case of a spherical infinitesimal cavity, the problems in [14] have an analytical solution, from which the following closed-form expression of A is obtained:

A ijk = 3(1-ν) 4μ(7-5ν) 5(δ ik δ j + δ jk δ i ) - 1+5ν 1+ν δ ij δ k [19]
Expression [18] of the topological derivative is implicit in that it relies upon solutions of boundary-value problems on the cavity-free reference body Ω. An explicit expression for T (x o ) is obtained in terms of the elastodynamic Green's tensor ûk (ξ, x, ω), i.e. the elastodynamic displacement field generated by a unit point force applied at x ∈ Ω along the k-th direction and satisfying the boundary conditions

ûk i (ξ, x, ω) = 0 (ξ ∈ S D , ξ = x), tk i (ξ, x, ω; n) = 0 (ξ ∈ S N , ξ = x)
In this case, the free and adjoint displacement fields are given by

u k (x) = S N p i (ξ)û k i (ξ, x, ω) dΓ ξ , u k (x) = S N p i (ξ)û k i (ξ, x, ω) dΓ ξ [20]
and the expression [18] becomes explicit as well. Such Green's tensors are explicitly known for simple unbounded media, but not in general for finite domains.

Numerical examples

Pressurized annular sphere

To validate the foregoing developments, the elastodynamic Neumann problem for a spherical shell (outer radius R, inner radius a true < R), shown in Fig. 1, is considered. The shell, Ω true , is subjected to a uniform time-harmonic pressure p acting over its external surface S. For this problem, the closed-form expression for topological derivative, T (r o ), where 0< r o <R denotes the radial coordinate of a sampling point inside the void-free, i.e. reference body Ω (see Fig. 1), can be obtained for the least squares cost function J with S obs = S and W = I 2 . To provide a basis for comparison, numerical values of the topological derivative T (r o ) are computed by means of a three-dimensional boundary element solution [START_REF] Pak R | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF] applied to the adjoint-field formula [18] with p = μ, ν = 0.3, R/a true = 3 and ωR ρ/μ = 3. In the numerical model, the inner and outer surfaces of the annular sphere were discretized using 486 and 150 eight-noded boundary elements, respectively, which provides at least 12 element lengths per shear wavelength for the excitation frequency chosen. The adjoint field technique for computing T (r o ) yields results that are in good agreement with the closed-form solution (see Fig. 1), the RMS error for the featured set of grid points being 1.5 10 -3 . It should be mentioned, however, that the accuracy of numerical estimates of T (r o ) was found to deteriorate for values of r o /R close to unity, primarily as a result of the near-singular nature of the integral representation for stresses at observation points close to the boundary.

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¡¢

Cavity embedded in a half-space

The configuration is as depicted in Fig. 2, with the 'true' spherical cavity, of diameter D = 0.4d, centered at (d, 0, 3d). In succession, the cavity is illuminated by 16 axial point sources acting on the surface of a semi-infinite solid. For each source location, the surface motion u obs is monitored at 25 sensors distributed over the square testing grid; here, this data is simulated using the BEM formulation of the forward problem. Four excitation frequencies ω ≡ ωd ρ/μ = 1, 2, 4, 8 have been considered.

For this testing configuration, the values of T (x o ) are computed over the horizontal surface S h = {ξ ∈ Ω| -5d < ξ 1 < 5d, -3d < ξ 2 < 3d, ξ 3 = 3d} passing through the centroid of the 'true' cavity and plotted in Fig. 3 for the above-defined set of frequencies. The computational grid is chosen so that the sampling points x o are spaced by 0.25d in both ξ 1 and ξ 2 directions. In the display, the red tones indicate negative values of T and thus possible cavity location; for comparison, the true cavity is outlined in white in each of the diagrams. The results clearly demonstrate the usefulness of the topological derivative as a computationally efficient tool for exposing the approximate cavity location, with "higher" frequencies (ω = 2, 4) providing in general better resolution. From the diagram for ω = 8 where λ s /D ≈ 1, however, it is also For completeness, the variation of T (x o ) across the vertical planar region S v = {ξ ∈ Ω|-5d < ξ 1 <5d, ξ 2 = 0, 0.25d < ξ 3 <6d} is given in Fig. 4. Similar to the earlier diagram, the sampling points x o are spaced by 0.25d in the ξ 1 and ξ 3 directions. The fact that both source and receiver points are limited to a single planar surface clearly leads to a diminished resolution compared to the previous result. The contour plots for ω = 2 and 4 exhibit greater accuracy than that for ω = 1, but are also plagued with local minima that are absent in the former diagram. The non-informative distribution With diagrams such as those in Figs. 3 and4, an algorithm for identifying plausible cavity locations could be thus devised on the basis of the non-zero distribution of

T (x o ) = T (x o )H(-C -T (x o )) [21]
where

Cavity embedded in a cube

To illustrate the utility of topological derivative as a preliminary tool for elasticwave sounding of finite bodies, this example deals with the delineation of a spherical In sequence, a virtual elastodynamic experiment is performed on each of the four vertical faces of the cube. With reference to the vertical testing face ξ 1 = -3d, the experiment parameters are chosen so that i) Γ obs = {ξ ∈ S N , |ξ 1 = -3d, -3d<ξ 2 < 3d, -3d < ξ 3 < 3d}, and ii) the cavity is illuminated in sequence by five localized (pyramid-shaped) distributions p q (ξ) (ξ ∈ Γ obs , q = 1, 2, . . . 5) of surface tractions, each of resultant P = μd 2 , centered at (-3d, 0, 0) and (-3d, ±1.8d, ±1.8d), respec- Figure 7 illustrates the distribution of topological derivative corresponding to excitation frequencies ω = 1 and ω = 2. In both displays, the distribution is plotted in three coordinate (cutting) planes containing the center of the true cavity. To provide a reference, intersection of each plane with the true cavity surface is outlined in white. One may observe that each distribution indeed points toward the true cavity through negative values of T (x o ). Unfortunately, the "low-frequency" map suffers from blurring, while an increased resolution of the "higher-frequency" map is diminished by the appearance of spurious minima. One should mention, however, that both excitation frequencies fall into the so-called resonance region [START_REF] Kress R | Inverse acoustic and electromagnetic scattering theory[END_REF] where the wave lengths are larger than the size of the scatterer: (λ s /d)| ω=1 = 2π and (λ s /d)| ω=2 = π, where d is the cavity diameter and λ s denotes the shear wave length. The particular effectiveness of long wavelengths for preliminary imaging by way of topological derivative is not surprising, since the assumption of an infinitesimal cavity, implicit to [9], is better conformed with by finite cavities that are 'small' relative to the probing wavelength.

ξ 1 /d ξ 1 /d ξ 2 /d ξ 3 /d c) ξ 2 /d ξ 2 /d ξ 3 /d a) b) ξ 1 /d ξ 1 /d ξ 2 /d ξ 3 /d c) ξ 2 /d ξ 3 /d

Summary

In this study, the concept of topological derivative that has its origins in elastostatics and shape optimization is extended to 3D elastic-wave imaging of solids. The proposed derivation follows a pattern which is generic, i.e. transposable to i) other numerical techniques (e.g. BIE, BEM, or FDM), ii) variety of physical problems described by linear field equations such as acoustics, heat transfer, and electromagnetics, and iii) other types of infinitely small objects, such as inhomogeneities or cracks. Numerical results suggest that the topological derivative is a useful tool for outlining hidden objects from surface measurements at a modest computational cost.

Figure 1 .

 1 Figure 1. Spherical shell: geometry and notation (top); topological derivative T (r o ) along a radial line (bottom)

Figure 2 .

 2 Figure 2. Cavity embedded in a half-space: definition and notation

Figure 3 .Figure 4 .

 34 Figure 3. Cavity embedded in a half-space: distribution of (μd) -1 T (x o ) in the ξ 3 = 3d (horizontal) plane

  C < 0 denotes a suitable threshold value and H(•) is the Heaviside step function. With such definition, it is also possible to combine the individual advantages of different probing wavelengths by employing the product of [21] at several frequencies. As an illustration of the latter approach, Fig. 5 plots the distribution of the product of T | ω=1 and T | ω=2 in the vertical plane, with C set to approximately 40% of the global minima of the respective distributions in Fig. 4. Despite the limited accuracy and multiple minima characterizing respectively the individual solutions for ω = 1 and ω = 2, the combined result stemming from [21] points to a single cavity with its centre and size closely approximating the true void configuration.

0Figure 5 .

 5 Figure 5. Distribution of (μd) -2 T| ω=1 × T| ω=2 in the ξ 2 = 0 plane

Figure 6 .

 6 Figure 6. Cavity embedded in a cube: BEM mesh

Figure 7 .

 7 Figure 7. Distribution of T (x o )/d in coordinate planes: a) ξ 3 = d, b), ξ 2 = 1.5d, and c) ξ 1 = d containing the center of the true cavity (top: ω = 1, bottom: ω = 2)

  In addition, Fig.8plots the distribution of the product T (x o )| ω=1 × T (x o )| ω=2 in three coordinate planes containing the center of the true cavity, where T (x o ) is defined by[21] with C = 0.5 inf x o T (x o ). Again, this distribution clearly points to a single cavity with its center and size closely approximating the true void configuration.

Figure 8 .

 8 Figure 8. Distribution of (1/d 2 ) T| ω=1 × T| ω=2 in coordinate planes: a) ξ 3 = d, b), ξ 2 = 1.5d, and c) ξ 1 = d containing the center of the true cavity

  ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤

	p
	R
	S
	£ £
	p
	S

Acknowledgements

This research is currently supported by a CNRS-NSF cooperation grant. In addition, special thanks are extended to École Polytechnique for providing a Visiting Scientist appointment to B. Guzina, and to MTS Systems Corporation for funding visits by M. Bonnet at the University of Minnesota as a MTS Visiting Professor.