Mark Asch
email: mark.asch@u-picardie.fr

Marion Darbas

Jean-Baptiste Duval

NUMERICAL SOLUTION OF AN INVERSE INITIAL BOUNDARY VALUE PROBLEM FOR THE WAVE EQUATION IN THE PRESENCE OF CONDUCTIVITY IMPERFECTIONS OF SMALL VOLUME

Keywords: Mathematics Subject Classification. 35R30, 35L05, 65M60 Wave equation, exact controllability, inverse problem, finite elements, Fourier inversion

come

Numerical solution of an inverse initial boundary value problem for the wave equation in the presence of conductivity imperfections of small volume

Introduction

The localization of small imperfections is of great importance since there are numerous practical applications, particularly in the fields of medical imaging and nondestructive testing of materials. Generally, when we seek to localize an imperfection contained in a bounded domain, we need to solve an inverse problem for retrieving the geometry of the imperfection.

The determination of conductivity profiles from knowledge of boundary measurements has received a great deal of attention (see, for example, [START_REF] Alves | Boundary integral formulae for the reconstruction of imperfections of small diameter in an elastic medium[END_REF][START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF][START_REF] Ammari | Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter[END_REF][START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF]). However, the reconstruction of imperfections within a dynamical (i.e. time-dependent) framework has not been widely investigated. The present paper is a first attempt to implement an effective numerical method to determine the location of small-sized, conductivity imperfections inside a homogeneous medium from dynamical measurements on a part of the boundary. In imaging, this is known as "limited-view data".

The inverse problem considered in this paper is more complicated from the mathematical point of view and more interesting from the point of view of applications than the one solved in [START_REF] Ammari | Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter[END_REF][START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF] because one often cannot get measurements for all t or on the whole boundary, and so one cannot, by taking a Fourier transform in the time variable, reduce our dynamic inverse problem to the Helmholtz equation considered in [START_REF] Ammari | Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter[END_REF][START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF]. Previous numerical investigations have concentrated on time independent equations with full data -see for example, [START_REF] Ammari | Direct elastic imaging of a small inclusion[END_REF][START_REF] Ammari | A MUSIC-type algorithm for detecting internal corrosion from electrostatic boundary measurements[END_REF][START_REF] Asch | Numerical Localizations of 3D Imperfections from an Asymptotic Formula for Perturbations in the Electric Fields[END_REF]. This article is based on three sources: the Hilbert Uniqueness Method (HUM) formulated by Lions -see [START_REF] Lions | Contrôlabilité exacte, Perturbations et Stabilisation de Systèmes Distribués, Tome 1, Contrôlabilité exacte[END_REF]; the numerical application of the above method to exact boundary control of the wave equation by Glowinski et al. -see [START_REF] Glowinski | A numerical approach to the exact controllability of the wave equation (I). Dirichlet controls: Description of the numerical methods[END_REF]; the theoretical results of the problem of detection for the wave equation by Ammari -see [START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF]. This paper is structured as follows. We begin, in Section 2, by defining some notation and then formulate the inverse initial boundary value problem. We also address the identification procedure by the Fourier method. In Section 3, we recall the numerical algorithm for the Hilbert Uniqueness Method of Lions [START_REF] Lions | Contrôlabilité exacte, Perturbations et Stabilisation de Systèmes Distribués, Tome 1, Contrôlabilité exacte[END_REF]. We present the numerical method used for the dynamic detection problem in Section 4. Then numerical results obtained from simulations are shown in Section 5. Finally, some conclusions and perspectives are reported in the last section.

Dynamic localization theory

In this part, we recall the different results contained in the article of Ammari (see [START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF] for the full details) on which the numerical algorithm and the subsequent simulations will be based.

Some notation and presentation of the inverse problem

Some notation

Let Ω be a bounded, smooth subdomain of R d , d = 2, 3, with for simplicity, a smooth boundary ∂Ω, with n denoting the outward unit normal to ∂Ω. We suppose that Ω contains a finite number m of imperfections, each of the form z j + αB j , where B j ⊂ R d is a bounded, smooth domain containing the origin. This gives a collection of imperfections of the form B α = ∪ m j=1 (z j + αB j ). The points z j ∈ Ω, j = 1, ..., m, that define the locations of the imperfections are assumed to satisfy two distance conditions:

|z j -z l | ≥ d 0 > 0 ∀ j = l, dist(z j , ∂Ω) ≥ d 0 > 0 ∀ j. (2.1)
We also assume that α > 0, the common order of magnitude of the diameters of the imperfections, is sufficiently small so that these are disjoint and their distance to R d \ Ω is larger than d 0 /2.

Presentation of the inverse problem

Let γ 0 denote the conductivity of the background medium which, for simplicity, we shall assume is constant.

Let γ j be the constant conductivity of the j-th imperfection, z j + αB j . We define the piecewise constant conductivity

γ α (x) = γ 0 , if x ∈ Ω\ Bα , γ j , if x ∈ z j + αB j .
Consider the initial boundary value problem for the wave equation, in the presence of the imperfections,

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂ 2 u α ∂t 2 -∇ • (γ α ∇u α ) = 0 in Ω × (0, T ), u α = f on ∂Ω × (0, T ), u α | t=0 = u 0 , ∂u α ∂t t=0 = u 1 in Ω. (2.2)
Define u to be the solution of the wave equation in the absence of any imperfections, satisfying

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂ 2 u ∂t 2 -∇ • (γ 0 ∇u) = 0 in Ω × (0, T ), u = f on ∂Ω × (0, T ), u| t=0 = u 0 , ∂u ∂t t=0 = u 1 in Ω. (2.3)
Here T > 0 is a final observation time, and the initial conditions ϕ, ψ ∈ C ∞ ( Ω) and the boundary condition

f ∈ C ∞ (0, T ; C ∞ (∂Ω)
) are subject to the compatibility conditions at the initial time,

∂ 2l t f | t=0 = (γ 0 ) l (Δ l u 0 )| ∂Ω and ∂ 2l+1 t f | t=0 = (γ 0 ) l (Δ l u 1 )| ∂Ω , l = 1, 2, . . .
These ensure that the initial boundary value problem (2.3) has a unique solution in C ∞ ([0, T ] × Ω); see [START_REF] Evans | Partial Differential Equations[END_REF].

The transmission problem for the wave equation (2.2) has a unique weak solution u α ∈ C 0 (0, T ; H 1 (Ω)) ∩ C 1 (0, T ; L 2 (Ω)); see, for example, [START_REF] Lions | Contrôlabilité exacte, Perturbations et Stabilisation de Systèmes Distribués, Tome 1, Contrôlabilité exacte[END_REF]. Moreover, it can be shown that ∂uα ∂n | ∂Ω belongs to L 2 (0, T ; L 2 (∂Ω)).

We employ standard L 2 -based Sobolev spaces to measure function regularity. The notation H s denotes those functions which, along with all their derivatives of order less than or equal to s, are in L 2 . In particular, H 1 0 denotes the closure of C ∞ 0 in the norm H 1 . Sobolev spaces with negative indices are defined by duality, using an L 2 -inner product. We shall only need one such space, namely, H -1 , which is defined as the dual of H 1 0 .

The inverse problem is: determine the centers of the imperfections from the knowledge of boundary data, generated numerically,

∂u α ∂n on Γ c × (0, T ),
where Γ c ⊆ ∂Ω is a part of the boundary that must satisfy certain geometrical constraints, and (0, T ) is the finite time-interval during which the measurements are recorded. It is also possible to compute certain properties of the shapes of the imperfections B α , though this is not treated here. For this purpose, Ammari, in [START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF], has developed an asymptotic method based on appropriate averaging of dynamic boundary measurements, using particular background solutions as weights. These particular solutions are constructed, when measurements are available on all or a part of the boundary of the domain, by the control method of Section 3.

The first fundamental step in the design of the reconstruction method is the derivation of an asymptotic formula for ∂uα ∂νj | ∂(zj+αBj ) + on the outer boundaries of the imperfections, in terms of the reference solution u, the location z j of the imperfection z j + αB j , and the geometry of B j . The next step is the use of this asymptotic formula to derive integral boundary formulae with a judicious choice of test functions and is based on the geometrical control method and the solution of Volterra-type integral equations. We will show that these boundary integral formulae form the basis of very effective numerical identification algorithms, as was predicted in [START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF]. We note that an analogous approach may be applied to the full (time-dependent) Maxwell equations with small imperfections of different electric permittivity or magnetic permeability (or both) -see [START_REF] Asch | Numerical Localizations of 3D Imperfections from an Asymptotic Formula for Perturbations in the Electric Fields[END_REF] for the time-harmonic case.

The asymptotic formula and the identification procedure

Following [START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF], we let β ∈ C ∞ 0 (Ω) be a cutoff function such that β(x) ≡ 1 in a subdomain Ω of Ω that contains the imperfections B α . For an arbitrary η ∈ R d , we assume that we are in possession of the boundary measurements of

∂u α ∂n on Γ c × (0, T ) for u 0 (x) = e iη•x , u 1 (x) = -i √ γ 0 |η| e iη•x and f (x, t) = e iη•x-i √ γ0|η|t .
(2.4) This particular choice of plane wave data u 0 , u 1 and f implies that the background solution u of the wave equation (2.3) is given explicitly by

u(x, t) = u η (x, t) = e iη•x-i √ γ0|η|t in Ω × (0, T ).
(2.5)

Suppose that T and the part Γ c of the boundary ∂Ω are such that they geometrically control Ω, which roughly means that every geometrical optics ray, starting at any point x ∈ Ω at time t = 0, hits Γ c before time T at a non-diffractive point -see [START_REF] Asch | Geometrical aspects of exact boundary controllability for the wave equation -A numerical study[END_REF][START_REF] Bardos | Sharp sufficient conditions for the observation, control and stabilization of waves from the boundary[END_REF]. Then from [START_REF] Lions | Contrôlabilité exacte, Perturbations et Stabilisation de Systèmes Distribués, Tome 1, Contrôlabilité exacte[END_REF], it follows that for any η ∈ R d , we can construct by the Hilbert uniqueness method (HUM), a unique g η ∈ H 1 0 (0, T ; L 2 (Γ)) in such a way that the unique weak

solution w η in C 0 (0, T ; L 2 (Ω)) ∩ C 1 (0, T ; H -1 (Ω)) of the wave equation ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ 2 w η ∂t 2 -∇ • (γ 0 ∇w η ) = 0 in Ω × (0, T ), w η = g η on Γ c × (0, T ), w η = 0 on ∂Ω\ Γc × (0, T ), w η | t=0 = β(x)e iη•x ∈H 1 0 (Ω), ∂w η ∂t t=0 = 0 in Ω, (2.6
)

satisfies w η (T ) = ∂ t w η (T ) = 0.
Next, for any η ∈ R d , let θ η be the solution to the Volterra equation of the second kind:

⎧ ⎪ ⎨ ⎪ ⎩ ∂θη(x,t) ∂t + ¡ T t e -i √ γ0|η|(s-t) θ η (x, s) -i √ γ 0 |η| ∂θη(x,t) ∂t ds = g η (x, t) for x ∈ Γ c , t ∈ (0, T ), θ η (x, 0) = 0 for x ∈ Γ c . (2.7)
The existence and uniqueness of this θ η in H 1 (0, T ; L 2 (Γ c )), for any η ∈ R d , can be established using the resolvent kernel. Since g η ∈ H 1 0 (0, T ; L 2 (Γ c )), the solution θ η belongs, in fact, to H 2 (0, T ; L 2 (Γ c )). We also note from differentiation of (2.7) with respect to t that θ η is the unique solution of the ODE

∂ t θ η -θ η = e i √ γ0|η|t ∂ t e -i √ γ0|η|t g η for x ∈ Γ c , t ∈ (0, T ), θ η (x, 0) = ∂ t θ η (x, T ) = 0 for x ∈ Γ c . (2.8)
It also immediately follows from this observation that θ η belongs to

H 2 (0, T ; L 2 (Γ)) since g η ∈ H 1 0 (0, T ; L 2 (Γ)).
To identify the locations and certain properties of the small imperfections B α , we average the boundary measurements ∂uα ∂n | Γc×(0,T ) , using the solution θ η to the Volterra equation (2.7) or, equivalently, the ODE (2.8) as a function of η. The following result holds.

Theorem 2.1 (Ammari).

Let η ∈ R d , d = 2, 3. Let u α be the unique solution in C 0 (0, T ; H 1 (Ω)) ∩ C 1 (0, T ; L 2 (Ω))
to the wave equation (2.2) with

u 0 (x) = u 0 η (x) = e iη•x , u 1 (x) = u 1 η (x) = -i √ γ 0 |η| e iη•x and f (x, t) = f η (x, t) = e iη•x-i √ γ0|η|t .
Suppose that Γ c and T geometrically control Ω [START_REF] Bardos | Sharp sufficient conditions for the observation, control and stabilization of waves from the boundary[END_REF]; then we have

¢ T 0 ¢ Γ θ η ∂u α ∂n - ∂u ∂n + ∂ t θ η ∂ t ∂u α ∂n - ∂u ∂n = - ¢ T 0 ¢ Γ e i √ γ0|η|t ∂ t e -i √ γ0|η|t g η ∂u α ∂n - ∂u ∂n = α d m j=1 γ 0 γ j -1 e 2iη•zj M j (η) • η -|η| 2 |B j | + o(α d ), (2.9)
where θ η is the unique solution to the ODE (2.8), with g η defined as the boundary control in (2.6), and M j is the polarization tensor of B j , defined by

(M j ) k,l = e k • ¢ ∂Bj ν j + γ 0 γ j -1 ∂Φ j ∂ν j (y) y • e l ds j (y) . (2.

10)

Here (e 1 , ..., e d ) is an orthonormal basis of R d , ν j is the outward normal of the j-th imperfection and Φ j is the solution of a Laplace equation in the j-th imperfection.

See [START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF] for the full formulation and details of the proof. This theorem is fundamental for our identification problems.

Remark 2.1. In this result, the control function is a technical device that is used in an adjoint formulation based on integration-by-parts. This control does not have a physical interpretation. Its role is to reduce the identification problem to a quantity that can be easily calculated on the boundary.

Remark 2.2. The choice of plane-wave initial data is only one of many possible choices. If we choose other initial functions, we indeed obtain different imaging algorithms (see [START_REF] Ammari | Transient wave imaging with limited-view data[END_REF]) but the complexity of the inverse problem remains the same, since we still have to compute the same boundary control.

We can now proceed to describe our identification procedure which is based on Theorem 2.1. We neglect the asymptotically small remainder in the asymptotic formula (2.1) and define Λ α (η) by

Λ α (η) = - ¢ T 0 ¢ Γ e i √ γ0|η|t ∂ t e -i √ γ0|η|t g η ∂u α ∂n - ∂u ∂n (2.11) ≈ α d m j=1 γ 0 γ j -1 e 2iη•zj M j (η) • η -|η| 2 |B j | .
The function Λ α (η) is computed in the following way. First, we construct the control g η in (2.6) for given η ∈ R d . Then from the boundary measurements ∂uα ∂n | Γ×(0,T ) , and the explicit knowledge of ∂u ∂n , we form the weighted boundary integrals that appear in the above expression of Λ α (η).

Recall that the function e 2iη•zj is exactly the Fourier transform (up to a multiplicative constant) of the Dirac function δ -2zj (a point mass located at -2z j ). From Theorem 2.1, it follows that the function Λ α (η) is (approximately) the Fourier transform of a linear combination of derivatives of point masses, or

Λα (x) ≈ α d m j=1 L j δ -2zj (x),
where L j is a second order, constant coefficient, differential operator whose coefficients depend on the polarization tensor M j defined by (2.10) (see [START_REF] Ammari | Polarization and Moment Tensors: With Applications to Inverse Problems and Effective Medium Theory[END_REF][START_REF] Capdebosq | A review of some recent work on impedance imaging for inhomogeneities of low volume fraction[END_REF] for its properties) and Λα (x) represents the inverse Fourier transform of Λ α (η).

The method of reconstruction we propose here consists, as in [START_REF] Ammari | Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter[END_REF], of sampling values of Λ α (η) at some discrete set of points and then calculating the corresponding discrete inverse Fourier transform. After a rescaling by - 1 2 , the support of this discrete inverse Fourier transform yields the location of the small imperfections B α . In other terms, once Λ α (η) is computed from dynamic boundary measurements on Γ c , we calculate its inverse Fourier transform. Then the asymptotic formula (2.1) in Theorem 2.1 asserts that this inverse Fourier transform is a distribution supported at the locations {z j } m j=1 .

Once the locations are known, we may calculate the polarization tensors (M j ) m j=1 by solving an appropriate linear system arising from (2.1). These polarization tensors give ideas on the orientation and relative size of the imperfections. From the leading term of Λ α (η) given by (2.1), we cannot reconstruct more details of the shapes of the domains B j . Higher order terms in the asymptotic expansion of Λ α (η), with respect to α, are needed to reconstruct the domains B j with high resolution. The properties of the polarization tensors are described in [START_REF] Ammari | Polarization and Moment Tensors: With Applications to Inverse Problems and Effective Medium Theory[END_REF][START_REF] Capdebosq | A review of some recent work on impedance imaging for inhomogeneities of low volume fraction[END_REF].

The theoretical number of data (sampling) points needed for an accurate discrete Fourier inversion of Λ α (η)

follows from Shannon's sampling theorem. We need (conservatively) order (h/δ) 2 sampled values of Λ α (η) to reconstruct, with resolution δ, a collection of imperfections that lie inside a square of side h. Numerical experiments in [START_REF] Ammari | Boundary integral formulas for the reconstruction of electromagnetic imperfections of small diameter[END_REF] for the two-dimensional (time-independent) inverse conductivity problem confirm this sampling and seem to suggest that the method is quite stable with respect to noise in measurements and errors in the different approximations.

The HUM approach to exact boundary controllability

In this section we recall very briefly the numerical method that will be employed in order to solve the exact controllability problem. Full details of the method may be found in [START_REF] Asch | Geometrical aspects of exact boundary controllability for the wave equation -A numerical study[END_REF][START_REF] Castro | Numerical approximation of the boundary control for the wave equation with mixed finite elements in a square[END_REF][START_REF] Glowinski | Ensuring well posedness by analogy; Stokes problem and boundary control for the wave equation[END_REF][START_REF] Glowinski | Exact and approximate controllability for distributed parameter systems[END_REF][START_REF] Glowinski | A numerical approach to the exact controllability of the wave equation (I). Dirichlet controls: Description of the numerical methods[END_REF].

Formulation of the problem

We consider the wave equation with control on a part of the boundary,

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂ 2 t -Δ u = 0 in Q = Ω × (0, T ), u| t=0 = u 0 , ∂ t u| t=0 = u 1 in Ω, u = g on Σ c = Γ c × (0, T ), 0 on Σ \ Σ c = Γ \ Γ c × (0, T ), (3.1)
where Ω is a bounded domain of R d , d = 2, 3, with boundary Γ, and Γ c ⊂ Γ is the part of the boundary where the control is applied.

The problem of the exact boundary controllability is then: "Given T, u 0 , u 1 , can we find a control g on Σ c

such that the solution of (3.1) satisfies

u(x, T ) = ∂ t u(x, T ) = 0 in Ω?"
The answer is positive if one takes T sufficiently large and one controls on a set large enough to encounter every ray of geometric optics (see [START_REF] Asch | Geometrical aspects of exact boundary controllability for the wave equation -A numerical study[END_REF][START_REF] Bardos | Sharp sufficient conditions for the observation, control and stabilization of waves from the boundary[END_REF]). A systematic and constructive method for computing such a control, g, is provided by the Hilbert uniqueness method (HUM) introduced by Lions [START_REF] Lions | Contrôlabilité exacte, Perturbations et Stabilisation de Systèmes Distribués, Tome 1, Contrôlabilité exacte[END_REF].

Description of the HUM

We describe briefly the HUM for the control of the wave equation (3.1) from a part of the boundary (see [START_REF] Lions | Contrôlabilité exacte, Perturbations et Stabilisation de Systèmes Distribués, Tome 1, Contrôlabilité exacte[END_REF]).

Let

E = H 1 0 (Ω) × L 2 (Ω), E = H -1 (Ω) × L 2 (Ω)
and define the operator

Λ : E -→ E
as follows:

(1) Take e = e 0 , e 1 ∈ E and solve from t = 0 to t = T

⎧ ⎪ ⎨ ⎪ ⎩ ∂ 2 t -Δ ϕ = 0 in Q = Ω × (0, T ), ϕ| t=0 = e 0 , ∂ t ϕ| t=0 = e 1 in Ω, ϕ(x, t) = 0 on Σ = Γ × (0, T ). (3.2)
(2) Then solve (backwards) from t = T to t = 0

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂ 2 t -Δ ψ = 0 in Q = Ω × (0, T ), ψ(x, T ) = 0, ∂ t ψ(x, T ) = 0 in Ω, ψ = ∂ϕ ∂n on Σ c = Γ c × (0, T ), 0 on Σ\ Σ c = Γ \ Γ c × (0, T ). (3.3)
(3) Finally, define the operator Υ as

Υe = {ψ t (x, 0), -ψ(x, 0)} ∈ E .
We have the original theorem of Lions [START_REF] Lions | Contrôlabilité exacte, Perturbations et Stabilisation de Systèmes Distribués, Tome 1, Contrôlabilité exacte[END_REF]:

Theorem 3.1 (Lions). Operator Υ is linear and continuous from E onto E ; moreover, if T is sufficiently large

(> T min = 2 x -x 0 L ∞ (Ω) ) and if Γ c is of the type Γ (x 0 ) = {x|x ∈ Γ, (x -x 0 ) • n x > 0}
where x 0 ∈ R d is an arbitrary point and n x is the outward normal to Γ at x, then Υ is an isomorphism from E onto E .

Application of the HUM to the wave equation

We now apply Theorem 3.1 to the control of the wave equation (3.1). Suppose that

u 0 ∈ L 2 (Ω) , u 1 ∈ H -1 (Ω)
are given. Then,

(1) take f = u 1 , -u 0 -i.e. we identify u with ψ;

(2) solve Υe = f to obtain e 0 , e 1 , the initial data for the ϕ wave equation (3.1);

(3) solve the ϕ wave equation (3.1) forwards in time using e 0 , e 1 as initial data;

(4) calculate the normal derivative of the solution ϕ and set g = ∂ϕ ∂n | Σ0 ;

(5) solve the ψ wave equation (3.3) backwards in time using g as the boundary data;

(6) finally, set u = ψ, then since ψ(x, T ) = 0, ψ t (x, T ) = 0 was imposed, g (the boundary control) gives the exact boundary controllability with

u(x, T ) = ∂ t u(x, T ) = 0, ∀x ∈ Ω.
We remark that the operator Υ is symmetric and E-elliptic. These properties imply that, in step (2), Υe = f can be solved by a conjugate gradient algorithm. We will describe the numerical approximation of the control of the system (3.1) in Section 4.3.

Numerical discretization of dynamic localization

In this section, we present the numerical discretization of our identification problem. First of all, we give the steps of the reconstruction algorithm. Secondly, we briefly describe the finite element method and the time-stepping scheme used for solving the different wave equations contained in this detection problem. Then, we discuss the numerical discretization of the BiGrid HUM method, seen in Section 3. Finally, we present the implementation of the Fourier inversion method and apply it to some tests on the asymptotic formula (2.1).

Summary of the identification procedure algorithm

Let us give the steps of the procedure used by our identification algorithm. First, we suppose a finite number of imperfections, z j + αB j for j = 1, ..., m, with conductivities γ j . Then, for each η in a discrete set D of values,

(1) Compute the solution u α of the wave equation (2.2) by a finite element method to simulate the boundary data ∂u α ∂n on Γ c × (0, T ) for the inverse problem.

(2) Compute the quantity ∂u ∂n on Γ c × (0, T ) which is explicitly known from (2.5).

(3) Calculate the control g η of (2.6) via the BiGrid HUM method.

(4) Form the quantity Λ α (η) from the left hand side of (2.11) with a suitable quadrature formula.

Finally, apply the inverse Fourier transform to Λ α (η) (η in the Fourier space D) to compute Λα (x) (x in the physical space Ω). This leads to the solution of the numerical inverse problem, namely the localization of the centers of the imperfections.

Time-space discretization for the wave equation

In this subsection, we present the numerical resolution of the wave equation: we combine finite elements for the space discretization and finite differences for the time discretization. This will then be used for solving the control problem (2.6) for g η and the inhomogeneous wave equation (2.2) for u α .

We consider the following, general wave equation with a Dirichlet boundary condition, defined on Ω × (0, T ):

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂ 2 u ∂t 2 -∇ • (γ∇u) = f in Ω × (0, T ), u = g on ∂Ω × (0, T ), u| t=0 = u 0 in Ω, ∂u ∂t | t=0 = u 1 in Ω, (4.1)
with f ∈ L 2 (L 2 (Ω); (0, T )), g ∈ H 1 0 (∂Ω), Ω ⊂ R d , d = 2, 3.
We assume that γ ∈ L ∞ (Ω), γ(x) > 0, ∀x ∈ Ω, and that u 0 , u 1 ∈ H 1 (Ω). Besides, we consider the function g as differentiable in time and satisfying the necessary compatibility conditions with u 0 and u 1 . In practice, in order to work in the homogeneous space H 1 0 (Ω), we use a standard extension technique on the function g (see for instance [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]).

For the finite element mesh, we consider shape-regular meshes, T h , that partition the domain Ω into a set of disjoint elements {K}, such that Ω = ∪ K∈T h K. The elements are rectangles (in 2D) or parallelepipeds (in 3D), aligned with the coordinate axes. The diameter of an element K is denoted by h K and the mesh size h is given by h = max K∈T h h K . We define the following sub-space of H 1 0 (Ω) on the mesh,

V h = v h ∈ C 0 ( Ω); v h = 0 on ∂Ω and v h | K ∈ Q 1 , ∀K ∈ T h ,
where Q 1 defines the space of polynomials of degree at most one in each of the d variables, thus spanned by and N J respectively. An element v h of V h can be written as

∀x ∈ Ω, v h (x) = NI i=1 v h (x i )ω i (x),
where ω i are the Q 1 basis functions associated to the interior nodes. Moreover, we have

∂ 2 t v h (x, t) = i∈I v i (t)ω i (x), with v i (t) := ∂ 2 t v h (x i , t), i ∈ I.
We denote by P h the L 2 -projection onto V h . We obtain the linear, second-order system of ordinary differential equations for the discretization of (4.1) by the finite element method

M U h (t) + KU h (t) = B h 0 < t ≤ T, U h (0) = U 0 h and U h (0) = V 0 h given, (4.2)
where

U h (t) = (u i (t)) i∈I , U h (t) = (u h (x i , t)) i∈I , U 0 h = P h u 0 and V 0 h = P h u 1 . The mass matrix, M, has coefficients M ij = ¢ Ω ω i (x)ω j (x)dx, 1 ≤ i, j ≤ N I , the stiffness matrix, K, K ij = ¢ Ω γ(x)∇ω i (x)∇ω j (x)dx, 1 ≤ i, j ≤ N I ,
and B h ∈ C NI , the right hand side vector,

B i = ¢ Ω f (x i , t) ω i (x)dx - j∈J K ij g(x j , t) - j∈J M ij ∂ 2 g ∂t 2 (x j , t).
To discretize (4.2) in time, we use the fully implicit, unconditionally stable, Newmark time-stepping scheme (see [START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF][START_REF] Wood | Practical time-stepping schemes[END_REF]).

The BiGrid HUM method

This part is devoted to the computation of the boundary control g η of the wave equation in (2.6). This corresponds to the crucial step (3) of the reconstruction algorithm. It is by now well-known that numerical algorithms for the boundary control for the wave equation (2.6) may diverge when the mesh-size tends to zero. These numerical instabilities are due to the high-frequency spurious oscillations. Among the several remedies proposed in the literature [START_REF] Castro | Boundary controllability of a linear semi-discrete 1-D wave equation derived from a mixed finite element method[END_REF][START_REF] Castro | Numerical approximation of the boundary control for the wave equation with mixed finite elements in a square[END_REF][START_REF] Infante | Boundary observability for the space discretization of the one-dimensional wave equation[END_REF][START_REF] Zuazua | Boundary observability for the finite-difference space semi-discretizations of the 2-D wave equation in the square[END_REF], we choose to use the preconditioned BiGrid algorithm introduced by Glowinski in [START_REF] Glowinski | Ensuring well posedness by analogy; Stokes problem and boundary control for the wave equation[END_REF]. We have faithfully followed [START_REF] Asch | Geometrical aspects of exact boundary controllability for the wave equation -A numerical study[END_REF][START_REF] Glowinski | A numerical approach to the exact controllability of the wave equation (I). Dirichlet controls: Description of the numerical methods[END_REF] in order to formulate our algorithm. The main numerical novelty here is the generalization of the BiGrid algorithm to the 3D case. Theoretical results on the convergence of a BiGrid algorithm for the control of the wave equation are very recent. We can cite the article of Ignat and Zuazua [START_REF] Ignat | Convergence of a two-grid method algorithm for the control of the wave equation[END_REF] that proves the convergence in the case of a finite-difference space semi-discretizations of the 2D wave equation in the square. The Fourier analysis techniques that they have used can be adapted to any space dimension and to full approximation schemes (finite differences or structured meshes of finite elements), and hence to our problem. At the moment, we observe such behaviour numerically. Our BiGrid algorithm is stable with respect to the mesh-size for regular meshes and smooth data. As an illustration, we report some convergence errors (discrete norms) in Tables 1 and2 relative to the control of the wave equation on the unit cube using as initial data

• Table 1: u 0 (x, y, z) = 10 sin(πx) sin(πy) sin(πz), u 1 (x, y, z) = 0,

• Table 2: u 0 (x, y, z) = e -64[(x-0.5) 2 +(y-0.5) 2 +(z-0.5) 2 ] , u 1 (x, y, z) = 0.

Table 1. Oscillatory data on the unit cube with a partial control.

Mesh size

h = 1 8 h = 1 16 h = 1 32 h = 1 64 Iterations 7 7 6 5 u 0 -u 0 h L 2 (Ω) u 0 L 2 (Ω)
0.0238 0.0207 0.0103 0.0051

u h (T ) L 2 (Ω) u h (0) L 2 (Ω)
0.0467 0.0214 0.0135 0.0087 Table 2. Smooth exponential initial data on the unit cube with a partial control.

Mesh size

h = 1 8 h = 1 16 h = 1 32 h = 1 64 u 0 -u 0 h L 2 (Ω) u 0 L 2 (Ω)
1.093 0.0632 0.0328 0.0109

u h (T ) L 2 (Ω) u h (0) L 2 (Ω)
0.8060 0.1120 0.0645 0.0353 (see [START_REF] Asch | Geometrical aspects of exact boundary controllability for the wave equation -A numerical study[END_REF] for details on the stopping criterion). We observe a linear in h decay of the L 2 -error as expected from our Q 1 finite-element approximation.

A partial control on Γ c = {0 ≤ x ≤ 1, y = 0, 0 ≤ z ≤ 1} ∪ {x = 0, 0 ≤ y ≤ 1, 0 ≤ z ≤ 1} ∪ {x = 1, 0 ≤ y ≤ 1, 0 ≤ z ≤ 1} is

Fourier method

The localization of the centers, z j , j = 1, . . . , m, requires that we numerically compute the inverse Fourier transform of the boundary integral, Λ α (η) -see Section 2.2. A theoretical framework for the number of points needed is provided by Shannon's Theorem.

Theorem 4.1 (Shannon). If the Fourier transform of a function h(t) is zero for all frequencies greater than a critical frequency, f c (the Nyquist frequency), then the continuous function h(t) can be uniquely determined from a knowledge of its sampled values by the formula

h(t) = T ∞ n=-∞ h(nT ) sin 2πf c (t -nT ) π(t -nT )
where T = 1/(2f c ) and (1/T ) is known as the Nyquist sampling rate.

From this theorem we conclude the two following facts:

(1) If the imperfections are known to be contained within a square/cube of side length 2M, we must sample Λ α (η) with a step size Δη = 1/(2M ).

(2) If we sample in a domain |η| < η max , then the resolution will be at most δ = 1/(2η max ).

Thus we require N δ = (2M/δ) d , d = 2, 3, sampled values of Λ α (η) in order to reconstruct, at a resolution δ, a set of imperfections contained inside the square/cube of size 2M. This implies that we need (4M η max ) d points.

We now, following the analysis of [START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF], take a closer look at the quantity Λ α (η). For clarity, let us assume that d = 2 and let us examine a simplified form by assuming that our data is of the form

Λ(η) = e 2iz•η
for a single imperfection, z = (z 1 , z 2 ). Now we sample η in the domain [-η max , η max ] × [-η max , η max ], uniformly with n 2 points and define the discretization step-size as

Δη = 2η max n + 1 •
We obtain the discrete approximation,

Λ h (η) = e 2iz1(-ηmax+(k-1)Δη)+2iz2(-ηmax+(l-1)Δη) , 1 ≤ k, l ≤ n.
Applying the inverse discrete Fourier transform (IDFT) to this quantity, we obtain

1 n 2 n k=1 n l=1
e 2iz1(-ηmax+(k-1)Δη)+2iz2(-ηmax+(l-1)Δη) e 2πi( (k-1)

n (r-1))+2πi( (l-1) n (s-1))
for all 1 ≤ r, s ≤ n. Let us look at the first coordinate,

Λ1 = 1 n n k=1 e 2iz1(-ηmax+(k-1)Δη) e 2πi( (k-1) n (r-1)
) .

An elementary calculation yields

Λ1 = 1 n sin (2η max z 1 ) sin π z1Δη π + r-1 n •
For n large (n → ∞), Λ1 → 0 unless the quantity in the denominator equals zero. This will occur when sin kπ = 0, k = 0, ±1, . . . which means that z1Δη π + r-1 n must be close to an integer. So, once the sign of z 1 is fixed, in order to approach only one integer value (and thus observe a single peak in the Fourier transform), the following condition must hold:

z 1 Δη π ≤ 1 2 since 0 ≤ r-1 n ≤ 1.
For large values of η max , we begin to observe instabilities due to the exponential terms in the asymptotic formula that become highly oscillatory. In order to control these oscillations, we impose the truncation η ≤ (η * , η * ) T in a discrete 2-norm and set the quantity Λ α (η) = 0 for all other values of η as was done in [START_REF] Asch | Numerical Localizations of 3D Imperfections from an Asymptotic Formula for Perturbations in the Electric Fields[END_REF][START_REF] Vogelius | Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter[END_REF].

Test based on the asymptotic formula (2.1)

We recall the asymptotic formula (2.1), neglecting the higher order term,

Λ α (η) = α d m j=1 γ 0 γ j -1 e 2iη•zj M j (η) • η -|η| 2 |B j | , d = 2, 3. (4.3)
Let us illustrate through a 3-D case the inherent precision and robustness of this formula. We denote by N e the number of Fourier points for each component. We report in Figure 2 the successful localization of three imperfections via the Shannon sampling approach without truncation by taking N e = 128, η max = 64. Extensive tests of the formula with different sampling strategies can be found in [START_REF] Duval | Identification dynamique de petites imperfections[END_REF].

Figure 2. Shannon sampling approach without truncation, using the asymptotic formula, z 1 = (0.51, 0.64, 0.32), z 2 = (0.66, 0.32, 0.47) and z 3 = (0.39, 0.49, 0.61).

Table 3. Discretization parameters for all numerical simulations.

Parameter 2-D 3-D N 16, 32 14 h 1/N 1/N Δt h/ √ 2 h/ √ 3 T 3 > 2 √ 2 4 > 2 √ 3 10 -6
10 -6

Numerical simulations

We present here some numerical results of our identification algorithm. In all examples, we take the domain Ω to be the unit square [0, 1] 2 or the unit cube [0, 1] 3 . We use Q 1 finite elements for the space discretization and an implicit, unconditionally stable Newmark scheme for the time discretization (see Sect. 4.2). We also suppose that when the control is active on the entire boundary we take a final time, T, according to Table 3. The spatial discretization, h, and the corresponding time step Δt are also given in Table 3, as well as the stopping criterion of the BiGrid HUM method, . The cutoff function, β, is equal to 1 on a subdomain Ω = (0.2, 0.8) d , d = 2, 3.

The conductivity γ 0 of the background medium is equal to 1. The Fourier parameters depend on the geometry considered (number of imperfections and location of each imperfection) and the resolution required. In all the figures, we plot contours of the inverse Fourier transform, Λα (x).

Numerical localization of imperfections in 2D

We now present results of the numerical localization procedure, following the steps laid out in Section 4.1.

The Shannon sampling approach with truncation

In Figure 3 we show results for localization using the Shannon approach with truncation of the η-values. In both cases, we fix α = 0.04, N e = 128 and η max = 33. The truncation is quite severe, η * = 11, but we obtain satisfactory results. We have also applied various windowing methods for the computation of the inverse Fourier transform. These simulations can be found in the thesis [START_REF] Duval | Identification dynamique de petites imperfections[END_REF].

Limited-view data

Since the HUM method permits us to control on a part of the boundary, we can simulate cases where numerical data are available only on a subset of ∂Ω. In this case, the causality of the wave equation will imply the need of a longer control time T and the eventuality of an increase in the number of iterations required for the HUM algorithm to converge. A detailed presentation can be found in [START_REF] Asch | Geometrical aspects of exact boundary controllability for the wave equation -A numerical study[END_REF]. In Figure 4 we report This result is very important, since it illustrates the potential of our algorithm for treating more realistic cases, where full measurements are rarely available.

Numerical localization of imperfections in 3D

We now turn to the three-dimensional case. We consider the Shannon sampling approach and assume that we have either full or partial numerical boundary data.

The Shannon sampling approach with truncation

In the first example, we have taken one imperfection with equal abscissa and ordinates. The imperfection is situated at (0.25, 0.25, 0.25) with a radius α = 0.05. The conductivity of the medium is equal to 1 and that of the imperfection is equal to 10. The number of data (sampling) points for the calculation of Λ α is N e × N e × N e = 16 × 16 × 16 and the space step h = 1/12. The Fourier inversion is performed with η max = 19

and a truncation at η * = 9. In Figure 5 we compare the numerical result with that obtained from the asymptotic formula, with and without truncation.

The second test case has three imperfections, centered at z 1 = (0.39, 0.63, 0.32), z 2 = (0.71, 0.27, 0.43), z 3 = (0.51, 0.47, 0.66) . The number of data (sampling) points for the calculation of Λ α is N e × N e × N e = 64 × 64 × 64, with η max = 40 and a truncation at η * = 9. The spatial discretization step is h = 1/14, the radius of the imperfections α = 0.02 and the conductivity γ j = 10. In Figure 6 we observe that all three imperfections are well located. A 3D view of the localized imperfections is shown in Figure 7.

Finally, we confirm the robustness of our method by a test with four imperfections centered at z 1 = (0.66, 0.32, 0.47), z 2 = (0.55, 0.71, 0.39), z 3 = (0.39, 0.63, 0.31), z 4 = (0.71, 0.42, 0.74) . We fix the parame-

ters N e = 64, η max = 40 and η * = 9. The radius of the imperfections is now α = 0.01 and the conductivity γ j = 10. In Figures 8 and9, we can observe that all four imperfections are well located. Figure 8. Localization in 3D of four imperfections, centered at z 1 = (0.66, 0.32, 0.47), z 2 = (0.55, 0.71, 0.39), z 3 = (0.39, 0.63, 0.31), z 4 = (0.71, 0.42, 0.74) . Views from the x-, y-and z-directions.

Limited-view data

We consider again the second test of the previous section, using now partial measurements on Γ

c = {0 ≤ x ≤ 1, y = 0, 0 ≤ z ≤ 1} ∪ {x = 0, 0 ≤ y ≤ 1, 0 ≤ z ≤ 1} ∪ {x = 1, 0 ≤ y ≤ 1, 0 ≤ z ≤ 1}.
In other words, measurements are assumed to be available on only three adjacent faces of the cube. In Figure 10 we notice that the imperfections are still well located. Relative errors for the localization of the centers are given in Table 4. These results illustrate both the accuracy and the robustness of the method. Note that the theoretical resolution is δ = 1/2η max = 0.0125.

Some remarks on convergence

To end this section, we discuss various aspects related to the convergence of our numerical identification algorithm. To fix ideas, we took a single imperfection, centered at z = [0.5, 0.5], in the domain [0, 1] × [0, 1]. We considered two meshes, 20 × 20 with h = 0.05 and 40 × 40 with h = 0.025. For each mesh we ran the detection algorithm for imperfections of radius α = 0.1, 0.05, 0.025, 0.0125. We define the error as e = zz c / z , where z c is the computed center. We first consider the convergence behaviour of the localization error as a function of the size of the imperfections -see Table 5 (left). Our computations reveal that as long as the mesh size is sufficiently small, we can accurately detect imperfections of any size. We start to lose precision as soon as the imperfection radius is of the order of the mesh size. However, the imperfection center is still accurately recovered. This can be explained by the robustness of the Fourier approach where we can increase the number of points in the Fourier domain and thus reduce the step-size while still maintaining a large enough bandwidth in order to avoid truncation effects. The computations show that the only critical parameters are η max , N e and η * . In other words, by adjusting these three Fourier parameters, we can obtain an equally good precision, e, for the two meshes and for all imperfection radii -see Table 5 (right) and Table 6. The convergence of our implementation of the HUM algorithm has already been discussed in Section 4.3.

Concerning the convergence of the complete identification algorithm with respect to the (space-time) mesh size, we observe that this is related to a complex interplay between the HUM algorithm, the FFT and the numerical integration used to evaluate the boundary measurements. In particular, we note that numerical FFT precision is highly dependent on the effects of aliasing, truncation and correct sampling in general (see Sect. 4.4

and [START_REF] Brigham | The fast Fourier transform and its applications[END_REF]). Thus we conclude that if the FFT and the numerical integration are supposed accurate enough, the overall convergence will be that of the HUM algorithm i.e. first order in space and time. • coarse grain, where we exploit the MPI commands for looping on chunks of different values of η;

• and fine grain, where the linear solvers are parallelized using the PETSc package [START_REF] Balay | PETSc Web page[END_REF].

In the first case, the speed-up should be linear in the number of processors. In the second case, we expect the speed-up to be superlinear. We present in Figure 11 a graph of CPU time versus the number of processors. For a sequence of 2D computations, we have taken h = 1/32 and N e = 16, thus yielding 1024 sample points in x-y space, 135 time steps and 256 sample points in Fourier space.

We remark that we obtain a perfectly linear speed-up in the number of processors.

Conclusion and perspectives

This numerical study has proven the feasibility and the robustness of the dynamic localization procedure proposed by Ammari in [START_REF] Ammari | An inverse initial boundary value problem for the wave equation in the presence of imperfections of small volume[END_REF] that is based on the wave equation. The two central parts of this procedure are the computation of a control function, via a geometric control method, and a Fourier inversion of a boundary integral. We have carefully implemented these two parts. The first using an algorithm based on the HUM (Hilbert Uniqueness Method) with a BiGrid, finite-element spatial approximation and a Newmark scheme in time, and the second using Shannon's theorem and windowing functions. Our implementation of the HUM also permits the study of limited-view data where the observations are available on only a part of the frontier of the object/domain that we are imaging.

The results obtained, in 2D and in 3D, are remarkably accurate. We can localize an arbitrary number of randomly placed imperfections with good precision. Particularly encouraging are the good results obtained in the limited-view cases that show that if we can construct the control function accurately enough, we can obtain the same localization precision with limited-view data as that which can be obtained from full boundary data.

The perspectives are numerous. First of all, we intend to extend the treatment to general geometries.

However, here we need to ensure the applicability of the BiGrid method which has been based here on a structured mesh. In order to overcome this, we could resort to uniformly controllable approaches as in [START_REF] Asch | Uniformly controllable schemes for the wave equation on the unit square[END_REF],

or spectral approaches [START_REF] Lebeau | Experimental study of the HUM control operator for linear waves[END_REF], that can both be used with unstructured meshes. A second perspective is to treat coupled photoacoustic [START_REF] Ammari | Mathematical Modelling in Photo-Acoustic Imaging[END_REF] and magnetoacoustic [START_REF] Ammari | Mathematical models and reconstruction methods in magnetoacoustic imaging[END_REF] approaches in the dynamic regime. There is also the possibility of extending the complete localization procedure to Maxwell's equations, for which the theory exists (see [START_REF] Ammari | Identification of small amplitude perturbations in the electromagnetic parameters from partial dynamic boundary measurements[END_REF]) but the numerical implementation of the geometric control still remains to be done.

Figure 1 .

 1 Figure 1. An example of a 2-D domain containing three imperfections, z 1 + αB 1 , z 2 + αB 2 , z 3 + αB 3 .

  {1, x, y, xy} when d = 2, and by {1, x, y, z, xy, xz, yz, xyz} when d = 3. We denote by I the set of interior nodes (the nodes not belonging to Γ = ∂Ω) and by J the set of boundary nodes. Their cardinalities are N I

  applied. In our calculations, we fix the time step Δt = h/ √ 3, T = 4 and = 10 -6

(a) z 1 =Figure 3 .

 13 Figure 3. Localization in 2D of three imperfections, with thresholding at η ≤ ||(11, 11)|| l 2 . The true centers are marked by the small squares.

Figure 4 .

 4 Figure 4. Localization in 2D with limited-view data.

Figure 5 .

 5 Figure 5. Localization in 3D of a single imperfection at z = (0.25, 0.25, 0.25): asymptotic formula (left), asymptotic formula with truncation (center), numerical localization (right)views of x-y plane at z = 0.25.

Figure 6 .

 6 Figure 6. Localization of three imperfections, centered at z 1 = (0.39, 0.63, 0.32), z 2 = (0.71, 0.27, 0.43), z 3 = (0.51, 0.47, 0.66) . Views from the x-, y-and z-directions.

Figure 7 .

 7 Figure 7. Localization of three imperfections centered at z 1 = (0.39, 0.63, 0.32), z 2 = (0.71, 0.27, 0.43), z 3 = (0.51, 0.47, 0.66) . 3D view.

Figure 9 .

 9 Figure 9. Localization of four imperfections centered at z 1 = (0.66, 0.32, 0.47), z 2 = (0.55, 0.71, 0.39), z 3 = (0.39, 0.63, 0.31), z 4 = (0.71, 0.42, 0.74) . 3D view.

Figure 10 .

 10 Figure 10. Localization in 3D with limited-view data of three imperfections, centered at z 1 = (0.39, 0.63, 0.32), z 2 = (0.71, 0.27, 0.43), z 3 = (0.51, 0.47, 0.66) . Views from the x-, y-and z-directions.

Remark 5 . 1 .

 51 In order to reduce the CPU time needed by our identification algorithm, we have used Message Passing Interface (MPI) commands for a parallel implementation in Fortran 90. Two levels of parallelism are possible:

Figure 11 .

 11 Figure 11. Performance curve for the localization of one imperfection with parallelization.

Table 4 .

 4 Localization of the centers of three imperfections. Relative localization errors (total control versus partial control).

	Theoretical centers	Total control Partial control
	z 1 = (0.39, 0.63, 0.32)	0.072	0.072
	z 2 = (0.71, 0.27, 0.43)	0.007	0.007
	z 3 = (0.51, 0.47, 0.66)	0.033	0.052

Table 5 .

 5 Convergence of the Fourier detection algorithm as a function of the imperfection size (on left) and of η max (on right) for the 20 × 20 mesh with h = 0.05.

	α	N e η max η *	e	α N e η max η *	e
	0.1 128 80	6 0.0149	0.1 64	40	6 0.0526
	0.05 128 80	6 0.0149	0.1 64	50	6 0.0301
	0.025 128 80	6 0.0149	0.1 64	60	6 0.0150

Table 6 .

 6 Convergence of the Fourier detection algorithm as a function of the imperfection size for the 40 × 40 mesh with h = 0.025.

	α	N e η max η *	e
	0.1	128 80	6 0.0149
	0.05 128 80	6 0.0149
	0.025 128 80	6 0.0149
	0.0125 128 80	6 0.0149

Acknowledgements. The authors would like to thank the referees whose highly relevant remarks greatly improved the final form of this paper.


Paulo B De Freitas
email: defreitaspb@gmail.com

Sandra M S F Freitas

Marcos Duarte

Mark L Latash

Vladimir M Zatsiorsky

Effects of joint immobilization on standing balance

Keywords: posture, quiet stance, inverted pendulum, joint constraint, postural control PsycINFO classification: 2300 Human Experimental Psychology Classification code: 2330 Motor Processes

We investigated the effect of joint immobilization on the postural sway during quiet standing. We hypothesized that center of pressure (COP), and rambling and trembling trajectories would be affected by joint immobilization. Ten young adults stood on a force plate during 60 s without and with immobilized joints (only knees constrained, CK; knees and hips, CH; and knees, hips, and trunk, CT), with their eyes opened (OE) or closed (CE). The root mean square deviation (RMS, the standard deviation from the mean) and mean speed of COP and the rambling and trembling trajectories in the anterior-posterior and medial-lateral directions were analyzed. Similar effects of vision were observed for both directions: larger amplitudes for all variables were observed in the CE condition. In the anterior-posterior direction, postural sway increased only when the knees, hips, and trunk were immobilized. For the medial-lateral direction, the RMS and the mean speed of the COP and the rambling and trembling displacements decreased after immobilization of knees and hips and knees, hips, and trunk. These findings indicate that the inverted pendulum model is unable to completely explain the processes involved in the control of the quiet upright stance in the anterior-posterior and medial-lateral directions.

Introduction

During upright quiet standing the body has often been modeled as a single inverted pendulum, pivoting about the ankle joints [START_REF] Gage | Kinematic and kinetic validity of the inverted pendulum model in quiet standing[END_REF][START_REF] Winter | Stiffness control of balance in quiet standing[END_REF]. The premise of the inverted pendulum model is that the body above the ankles behaves as a rigid structure. The use of the inverted pendulum model to represent balance maintenance in quiet standing has been supported by kinematic, kinetic [START_REF] Gage | Kinematic and kinetic validity of the inverted pendulum model in quiet standing[END_REF][START_REF] Karlsson | The ankle strategy for postural control: A comparison between a model-based and a marker-based method[END_REF], and EMG [START_REF] Gatev | Feedforward ankle strategy of balance during quiet stance in adults[END_REF] measures in both the sagittal and frontal plane [START_REF] Winter | Unified theory regarding A/P and M/L balance in quiet stance[END_REF]. However, most of the kinematic studies have shown small angular displacement at joints above the ankle [START_REF] Gage | Kinematic and kinetic validity of the inverted pendulum model in quiet standing[END_REF][START_REF] Gatev | Feedforward ankle strategy of balance during quiet stance in adults[END_REF]. For example, Gage and collaborators [START_REF] Gage | Kinematic and kinetic validity of the inverted pendulum model in quiet standing[END_REF] reported that knee joint movements were larger as compared to ankle joint displacement in quiet standing but they still considered the inverted pendulum model a good representation of quiet standing in humans. Conversely, it has also been claimed that small angular displacements at joints above the ankles play an important role in the control of quiet upright stance and their contribution should not be ignored (Aramaki, Nozaki, Masani, Sato, Nakazawa, & Yano, 2001;[START_REF] Bardy | Postural coordination modes considered as emergent phenomena[END_REF][START_REF] Creath | A unified view of quiet and perturbed stance: simultaneous co-existing excitable modes[END_REF][START_REF] Day | Effect of vision and stance width on human-body motion when standing: Implications for afferent control of lateral sway[END_REF][START_REF] Fitzpatrick | Stable human standing with lower-limb muscle afferents providing the only sensory input[END_REF][START_REF] Krishnamoorthy | Joint coordination during quiet stance: Effects of vision[END_REF][START_REF] Nicholas | A link-segment model of upright human posture for analysis of head-trunk coordination[END_REF].

Recently, Hsu and collaborators [START_REF] Hsu | Control and estimation of posture during quiet stance depends on multijoint coordination[END_REF] and Krishnamoorthy and colleagues [START_REF] Krishnamoorthy | Joint coordination during quiet stance: Effects of vision[END_REF] reported that stabilization of a reference position (assumed to be described by the position of the center of mass, COM) was achieved by coordination of several joints along the human ''pendulum'' when participants stood quietly on a normal and on a narrow base of support. Similar findings were observed when sudden perturbations were applied to the participants in quiet stance (Alexandrov, Frolov, Horak, Carlson-Kuhta, & Park, 2005), when supra-postural tasks were performed in quiet standing [START_REF] Bardy | Postural coordination modes considered as emergent phenomena[END_REF][START_REF] Bardy | Dynamics of human postural transitions[END_REF][START_REF] Bardy | On perturbation and pattern coexistence in postural coordination dynamics[END_REF], during fast-forward voluntary trunk bending movements (Alexandrov, Frolov, & Massion, 2001), and during whole-body voluntary movements [START_REF] Freitas | Two kinematic synergies in voluntary whole-body movements during standing[END_REF]. In addition, [START_REF] Kuo | Human standing posture: Multijoint movement strategies based on biomechanical constraints[END_REF] found that the "ankle strategy" involving rotation about the ankles only [START_REF] Nashner | Analysis of stance posture in humans[END_REF] requires more muscle activation than the "hip strategy" (two-segment inverted pendulum pivoting about the hip and ankle joints) for a given amount of horizontal acceleration. Hence, the hip strategy is more effective at controlling the COM position with minimal muscle activation when comparing to the "ankle strategy".

Another approach to investigate the adequacy of the single pendulum model to describe the control of the upright stance and, concurrently, the role played by the joints above the ankles in this control is to prevent the movement of those joints while maintaining a quiet stance. To date, only two studies investigated the effect of preventing movements of the joints above the ankles on postural control (Aramaki et al., 2001;[START_REF] Fitzpatrick | Stable human standing with lower-limb muscle afferents providing the only sensory input[END_REF], and only Fitzpatrick and colleagues (1994) measured the amount of postural sway in the anterior-posterior direction after joint immobilization, whereas the effect of joint immobilization on medial-lateral direction was not assessed.

They observed that when a rigid splint prevented the movement of all joints above the ankles, the magnitude of the ankle angular sway increased significantly.

The most common measure of postural sway in quiet standing is the center of pressure (COP) migration derived from force plate data. COP is the point of application of the resultant of vertical forces acting on the surface of support [START_REF] Duarte | Effects of body lean and visual information on the equilibrium maintenance during stance[END_REF]. Several methods have been proposed to decompose the COP displacement, commonly called stabilogram, with the intention to understand the mechanisms of postural control during quiet standing (Baratto, Morasso, Re, & Spada, 2002;[START_REF] Bottaro | Body sway during quiet standing: Is it the residual chattering of an intermittent stabilization process?[END_REF][START_REF] Collins | Open-loop and closed-loop control of posture: A random-walk analysis of center-of-pressure trajectories[END_REF][START_REF] Dijkstra | A gentle introduction to the dynamic set-point model of human postural control during perturbed stance[END_REF][START_REF] Frank | Multivariate Ornstein-Uhlenbeck processes with mean-field dependent coefficients: Application to postural sway[END_REF][START_REF] Winter | Stiffness control of balance in quiet standing[END_REF][START_REF] Zatsiorsky | Instant equilibrium point and its migration in standing tasks: rambling and trembling components of the stabilogram[END_REF][START_REF] Dijkstra | A gentle introduction to the dynamic set-point model of human postural control during perturbed stance[END_REF]. [START_REF] Zatsiorsky | Instant equilibrium point and its migration in standing tasks: rambling and trembling components of the stabilogram[END_REF] proposed a method to decompose the COP trajectory based on the idea that the body sways due to two reasons: (a) migration of the reference point with respect to which the equilibrium is instantly maintained (rambling), and (b) oscillation of the body with respect to the moving equilibrium reference (trembling). They found that the trembling migration correlates highly (and negatively) with the horizontal ground reaction force while such a correlation does not exist for the rambling. These data were interpreted as evidence that the trembling represents a deviation from the instant equilibrium position corrected via the horizontal ground reaction forces inducing the center of mass acceleration in the direction opposite to the deviation while the central controller does not correct the COP migrations associated with the rambling. According to the rambling and trembling hypothesis, the rambling and trembling reflect two distinct processes contributing to the control of body sway. The rambling trajectory may be related to supraspinal processes that are involved in the control of the migration of the instant reference point. The trembling may be a result of both, the action of spinal reflexes and changes in the intrinsic mechanical properties of the muscles and joints. In particular, the authors [START_REF] Zatsiorsky | Instant equilibrium point and its migration in standing tasks: rambling and trembling components of the stabilogram[END_REF][START_REF] Dijkstra | A gentle introduction to the dynamic set-point model of human postural control during perturbed stance[END_REF] suggested that the trembling was due to the moment at the ankle joints during quiet standing.

The analysis of the trajectories of the COP and its components during quiet standing while the joints above the ankles are immobilized could provide additional information regarding the effectiveness of the single pendulum model and also the role played by the above-the-ankle joints in the control of psoture. Therefore, the purpose of this study was to examine the contribution of the main joints (ankle, knee, hip, and joints of the trunk) to the control of quiet upright stance. Specifically, the current study investigated the COP, and the rambling and trembling behaviors when the knees, hips, and lower and medium portion of the trunk were constrained and the body balance had to be controlled exclusively by the non-constrained joints, specifically the ankle joints. We hypothesized that if the body behaves as a single inverted pendulum, the joint constraints will have a minor or no effect on the COP, and the rambling and trembling variables, whereas significant changes in the COP, the rambling, and, principally, in the trembling would indicate that the inverted pendulum model cannot accurately explain the processes involved in the control of quiet upright stance.

Because the hip abductors/adductors, in addition to the ankle invertors/evertors, have been considered responsible by controlling the medial-lateral movements [START_REF] Winter | Unified theory regarding A/P and M/L balance in quiet stance[END_REF], it was expected that the hip joint constraints would affect body sway in the frontal plane as well. Specifically, we predicted an increase in postural sway in medial-lateral direction mainly when hip and trunk were immobilized. In addition, this study investigated the combined effect of immobilization and vision on the sway in both directions (sagittal and frontal), while the participants stood with their eyes open or closed.

Method

Participants

Ten healthy adults, 5 males and 5 females, without any known neurological and musculoskeletal disorders, or balance problems, participated in this study. The participants' mean (±1 SD) age, body mass, and height were 28.3 (±4.6) years, 66.7 (±4.6) kg, and 169 (± 9) cm, respectively. Prior to testing the participants signed the informed consent form in accordance with the policies established by the Office for Regulatory Compliance of The Pennsylvania State University.

Procedures

The participants were instructed to stand upright as still as possible in bipedal stance on a 40×80 cm force platform (model 4080S, Bertec, Worthington, OH, USA). They stood barefoot with their feet positioned side-by-side at a comfortable width (about shoulder width), and arms placed in a relaxed position at the body's sides. Also, they were instructed to maintain the head erect and to look at a circular target (radius = 2.5 cm) positioned approximately 120 cm in front of them at eye level.

During the experiment the participants were submitted to four joint immobilization conditions: no constraint (NC), constraint of the knees (CK), constraint of the hips and knees (CH), and constraint of the trunk, hips, and knees (CT) and two visual conditions, open and closed eyes (OE and CE, respectively). Fig. 1 depicts the different levels of immobilization. The joint constraints were applied for each limb separately.

INSERT FIGURE 1 HERE

The joints were tightly immobilized by splints, casts, and duct tape. The knees constraints started from the superior part of the shank and ended at the inferior part of the thigh, with each knee immobilized separately. Four splints, one on the anterior, one on the posterior, one on the medial, and one on the lateral side, were placed vertically around each knee. After the splints were positioned, cast tapes were placed to fix the splints and finally duct tape was placed around each knee to increase the resistance and prevent any knee's flexion-extension movement. The hip was immobilized from of the proximal part of the thigh to of the inferior part of the trunk at the navel level. In the same way as the knee constraints, splints were placed around the thighs, pelvis and the inferior portion of the trunk; fixed by cast tapes and reinforced by duct tape. Finally, the trunk was immobilized from the navel line to the sternum xiphoid process, also using splints, casts, and duct tape. The immobilizations assured that joint movements around knees, hips, and trunk were completely restricted. During each trial the participants stood quietly for 60 s. Five participants started the experiment with the trunk, hips, and knee immobilized (CTOE and CTCE condition) and had the immobilization hardware removed in a top-to-down sequence (CT -NC). The other five participants started without any constraint; the immobilization was applied in a down-to-top sequence (NC -CT). This procedure was used to avoid fatigue effects in body sway.

Data acquisition and analysis

Data were recorded using a PC with a 12-bit A/D board (model AT-MIO-64E-3, National Instruments Corporation, Dallas, TX, USA) controlled by a special code written in LabView software (LabView 6.1; National Instruments, Dallas, TX, USA). The force platform used during the experiment provided information about the three forces and three moment components. The signals were used to calculate the center of pressure (COP) displacement in the anterior-posterior (COP AP ) and medial-lateral direction (COP ML ). The data were recorded with a sampling frequency of 100 Hz. The data analysis was performed using MATLAB software (MATLAB 6.5, The MathWorks, Natick, MA, USA). Data from horizontal forces (Fy and Fx) and COP (COP AP and COP ML ) were low-pass filtered with a zero-lag, fourth-order Butterworth filter with a cut-off frequency of 8 Hz.

Rambling and trembling were determined using the procedures described previously by [START_REF] Zatsiorsky | Instant equilibrium point and its migration in standing tasks: rambling and trembling components of the stabilogram[END_REF]. A brief description of the analyses is provided below. At the instances when the horizontal force (F hor ) is zero (F hor = 0) the body is instantaneously in an equilibrium state [START_REF] Zatsiorsky | An algorithm for determining gravity line location from posturographic recordings[END_REF]. The instances when F hor changed its sign were determined using local linear interpolation of the F hor time-series. Then, the instant equilibrium points (IEP, i.e., the COP locations at the instances when the horizontal forces were zero) were identified. After this procedure, the IEPs were used to compute the rambling and trembling trajectories. The IEPs and the rambling and trembling trajectories were determined separately from COP AP and COP ML time-series.

The rambling trajectory was established by interpolating the discrete IEPs with a cubic spline function. According to the rambling-trembling hypothesis, the rambling represents migration of the reference point on the supporting surface with respect to which the body equilibrium is maintained. The trembling represents the deviation of the COP from the equilibrium position and it was determined by subtracting the rambling trajectory from the COP trajectory (for more details see [START_REF] Zatsiorsky | Instant equilibrium point and its migration in standing tasks: rambling and trembling components of the stabilogram[END_REF][START_REF] Dijkstra | A gentle introduction to the dynamic set-point model of human postural control during perturbed stance[END_REF]. Fig. 2 depicts an example of the F hor trajectory, the IEPs, COP, and the rambling and trembling for a single trial in NCOE condition (anterior-posterior direction, 20 s).

INSERT FIGURE 2 HERE

The dependent variables analyzed in this study were root mean square deviation (RMS, the standard deviation from the mean) and mean speed of the COP, and the rambling and trembling displacement in the anterior-posterior (AP) and medial-lateral (ML) directions. Before calculating these variables, the very low frequencies were removed from the COP, and the rambling and trembling timeseries were linearly detrended by subtracting their own best fit line from their respective time-series. The RMS in both directions was determined by calculating the standard deviation of each signal, i.e., the oscillations of the signals around their mean values. The mean speed was determined by dividing the total distance along the signal trajectory by the total recording time.

Statistical analysis

Before performing the statistical analyses, we normalized, for each individual participant, the values of RMS and mean speed of the COP, and the rambling and trembling as percentage of their respective values measured in the NCOE condition (i.e., NCOE = 0%). For each variable, we subtracted the value of each one of the immobilization conditions by the value of the NCOE condition, divided the result by the value of the NCOE condition, and then multiplied by 100 (e.g., ((RMS CTOE -RMS NCOE )/ RMS NCOE )*100)). The statistical analyses were performed in SPSS (version 10.1, Chicago, IL) using the normalized values averaged over all trials for each participant for all dependent variables (RMS and mean speed of the COP, and the rambling and trembling time-series). To examine the effects of the joint immobilization (NC, CK, CH, CT) and vision (OE and CE) on each dependent variable, a total of four repeated measures two-way MANOVAs, two for AP and two for ML direction, were performed. When necessary, univariate analysis and post hoc tests with Bonferroni corrections were carried out. The significance level was set at .05. However, marginal effects (p < .1) were also considered in the analysis.

Results

All participants were able to stand with their joints immobilized. However, changes in postural sway due to joint immobilization were observed. Fig. 3 shows a stabilogram with a typical COP trajectory in each joint immobilization and visual condition. As can be seen, the COP displacements in the anterior-posterior (COP AP ) direction and in the medial-lateral (COP ML ) direction were influenced differently by the joint immobilization. These effects were observed for both visual conditions. Detailed statistical analyses of each variable of postural sway are presented below for each direction separately.

Table 1 presents the non-normalized means and standard errors values in both directions.

INSERT FIGURE 3 HERE INSERT TABLE 1 HERE

Anterior-posterior direction

The results showed that postural sway in anterior-posterior (AP) direction was altered after immobilization of knee, hips, and lower and middle trunk joints. Fig. 4

Medial-lateral direction

With respect to postural sway in medial-lateral (ML) direction, the results showed that when the participants had their knee and hip joints and their knees, hips, and trunk immobilized their postural sway in ML direction was altered.

Discussion

As expected, the results of the study, as summarized in Table 2, showed that postural sway was affected by the absence of visual information as well as joint immobilization. In addition, the presence or absence of visual information affected postural sway at each level of joint immobilization in a similar fashion. Overall, independent of COP direction (AP or ML), the rambling and trembling displacements were larger and quicker when participants kept their eyes closed as compared to eyes open. The results further showed that joint immobilization affects postural sway differently in the AP and ML directions.

In the AP direction, COP and trembling displacements and speeds were larger than in the control condition (NC) only when the knees, hips, and trunk were immobilized (CT). In contrast, in the ML direction, the COP and the rambling and trembling displacements and speeds dramatically reduced after immobilization of the hip joints (CH and CT conditions). Due to the distinctive effect of joint immobilization in AP and ML directions, the results for each direction will be discussed separately.

INSERT TABLE 2 HERE Anterior-posterior direction

The available literature data on the mechanics of the body sway in the anterior-posterior direction are controversial. On the one hand, the results of some studies have indicated that the body sway is controlled mainly at the ankle joints without any significant participation of other joints or muscle groups [START_REF] Gage | Kinematic and kinetic validity of the inverted pendulum model in quiet standing[END_REF][START_REF] Gatev | Feedforward ankle strategy of balance during quiet stance in adults[END_REF][START_REF] Karlsson | The ankle strategy for postural control: A comparison between a model-based and a marker-based method[END_REF][START_REF] Winter | Stiffness control of balance in quiet standing[END_REF]. By considering the body above the ankles as a rigid structure, such studies have proposed the inverted pendulum model for describing postural control during quiet standing. On the other hand, other studies have claimed that body sway is controlled by compensatory movements from several joints and not a single joint [START_REF] Krishnamoorthy | Joint coordination during quiet stance: Effects of vision[END_REF][START_REF] Peterka | Simplifying the complexities of maintaining balance[END_REF]. According to this view, all joints co-vary, that is, act as a synergy, trying to control the COM displacement (Alexandrov, Frolov, & Massion, 2001;[START_REF] Creath | A unified view of quiet and perturbed stance: simultaneous co-existing excitable modes[END_REF][START_REF] Fitzpatrick | Stable human standing with lower-limb muscle afferents providing the only sensory input[END_REF][START_REF] Freitas | Two kinematic synergies in voluntary whole-body movements during standing[END_REF], and minimize the overall joint variance [START_REF] Krishnamoorthy | Joint coordination during quiet stance: Effects of vision[END_REF], which consequently could reduce the muscle activity in a specific joint (for example, the ankle joint) [START_REF] Kuo | Human standing posture: Multijoint movement strategies based on biomechanical constraints[END_REF]. Thus, if this is the case, it could be expected an increase in postural sway when one or more joints were immobilized.

The results of the present study revealed that only after immobilization of the knees, hips, and trunk a small but significant increase in postural sway in AP direction could be identified. These results lead us to refute the hypothesis that postural stability in AP direction can be achieved without significant motion of joints above the ankles and to question the idea that the inverted pendulum model can adequately describe the upright posture during unperturbed quiet standing in terms of its whole body sway. These results are in accordance with the results of Fitzpatrick and colleagues (1994) who reported a 52% increase in body sway (i.e., ankle joint sway) when participants were splinted to limit the sway to the ankles than when the participants stood freely. These authors suggested that, when all body segments are free to move, there is an increase in body stability and that the movements around the ankle joints are not the most determinant factor in body stability.

It is important to mention that in the current study, unlike that of Fitzpatrick and colleagues (1994), we did not observe any change in RMS of the postural sway variable (i.e., COP). Whereas Fitzpatrick and colleagues (1994) assessed the RMS of the ankle angle sway (in radians) we assessed the RMS of the COP displacement (cm) in AP direction. However, the different variables have to be linked in the condition when all joints above the ankle are constrained. Therefore, a possible explanation for this discrepancy in results is that Fitzpatrick and colleagues (1994) immobilized the upper trunk, arms and head while we kept them unconstrained. Although the upper trunk, arms, and head only make a small contribution to the global COM/COP position, the obviation of head-trunk relative motion rather than lower limbs joint immobilization might have contributed to the increase in body sway. However, this is something yet to be investigated. Clearly, such lack of effect of joint mobilization is only possible in a quiet standing task, where the task offers no or little challenge to healthy individuals. It is well known that in a situation where the individual is perturbed during standing, stereotypical movements across all joints are usualy observed [START_REF] Nashner | Organization of rapid responses to postural and locomotorlike perturbations of standing man[END_REF].

The observed effect of immobilization on COP speed still leads us to believe that increased body stability is achieved when more joints are free to move. One explanation is that joints' immobilization would increase the level of activation of muscle groups located around the ankle joints which, consequently, could increase the stiffness around ankle joint and prevent increase in postural sway. The increase in RMS and mean speed of the trembling trajectories supports the idea that the trembling may be a result of the action of spinal reflexes and changes in the intrinsic mechanical properties of the muscles. However, we were not able to measure muscle activity and, then, more studies are necessary to support this hypothesis.

Medial-lateral direction

One of the goals of this study was to identify the effects of joint immobilization on medial-lateral (ML) direction. The results surprisingly revealed a reduction of RMS and speed of the COP and the rambling and trembling displacements in ML direction when knee and hips and knee, hips and lower and middle trunk joints were immobilized. According to Winter and collaborators [START_REF] Winter | Medial-lateral and anterior-posterior motor-responses associated with center of pressure changes in quiet standing[END_REF] in the ML direction the magnitude of body sway is constrained at the hip joints as a result of the activation of hip adductors and abductors. Therefore, one would expect that the immobilization of the knee joint would not affect the postural sway in the ML direction because most muscle groups around the knee are related to motions in the anterior-posterior direction. On the other hand, the immobilization of knee and hips and knee, hips and trunk would cause an increase in postural sway in ML direction, due the fact that the ankle invertors and evertors alone would not be able to restrict sufficiently the COP displacement. In contrast to our expectations, we observed the opposite effect, i.e., a decrease in COP and the rambling and trembling RMS and speed during the immobilization of knee and hips and knee, hips, and trunk. Carpenter and collaborators [START_REF] Carpenter | Surface height effects on postural control: A hypothesis for a stiffness strategy for stance[END_REF][START_REF] Carpenter | The influence of postural threat on the control of upright stance[END_REF] verified that when people maintained their upright stance in a threatening situation they decreased displacement amplitude and increased mean frequency due to augmentation of the stiffness caused by increased muscle co-contraction. However, our result can not be explained by the increased joint stiffness at the ankle joints, for instance by co-contraction of the ankle joint muscles, and, consequently, reduced postural sway due the fact that postural sway speed also decreased.

A possible explanation for the decrease in sway in the ML direction after joint immobilization would be that, with constraint hip movements, the individuals were less likely to use the mentioned load/unload mechanism resulting in lower COP displacement in the ML direction. This would be possible only if the displacement observed in the ML direction was not due to internally generated perturbations itself but due to self-motion and that we are also able to control our quiet standing posture in the ML direction without the use of the load/unload mechanism. The idea of self-motion to explain part of the body sway during quiet standing has indeed been proposed before as a possible exploratory mechanism [START_REF] Roebuck | Neuromuscular responses to disturbance of balance in children with prenatal exposure to alcohol[END_REF]. In addition, it might be possible that the load/unload mechanism using the hip joints is not an exclusive requirement for posture control in the ML direction during quiet standing, but only one of the strategies repertoire [e.g., the ankle eversion/inversion muscles can also participate in the control of standing [START_REF] Winter | Unified theory regarding A/P and M/L balance in quiet stance[END_REF][START_REF] Winter | Medial-lateral and anterior-posterior motor-responses associated with center of pressure changes in quiet standing[END_REF]. The load/unload mechanism using the hip joints has been shown to be essential for posture control in more challenging tasks involving weight transfers, such as gait initiation, walking, unconstrained standing, among others [START_REF] Freitas | Two kinematic synergies in voluntary whole-body movements during standing[END_REF]. In this last view, clearly less sway would not mean better control, only that a different mechanism of control is being used during quiet standing.

Another related interpretation for the contrasting effects of joint immobilization on the sway in the AP and ML directions allow for the following hypothesis. The ankle, knee, and hip joints show higher kinematic mobility in the AP direction resulting in a kinematically redundant system able to stabilize the center of mass location by co-varied changes in these joint angles [START_REF] Freitas | Two kinematic synergies in voluntary whole-body movements during standing[END_REF]. Motor redundancy may be viewed as a useful feature that allows, in particular, to deal with secondary tasks and external perturbations (Latash, Scholz, & Schöner, 2007;[START_REF] Zhang | What do synergies do? Effects of secondary constraints on multidigit synergies in accurate force-production tasks[END_REF]. Joint fixation may be viewed as an example of such a perturbation.

Apparently, the kinematic redundancy of the postural chain allowed to keep the AP sway characteristics unchanged under joint immobilization. In contrast, there is limited joint mobility in the ML direction and, as such, the system may be viewed as less (or even non-) redundant. For such a system, any additional constraint may lead to an impossibility of using the earlier mode of control. Note that, although our manipulations were designed to lead to joint immobilization primarily in the AP direction, they likely limited joint mobility in the ML direction as well, partly by limiting muscle length changes for postural muscles with the lines of action in-between the AP and ML axes. If we assume that sway magnitude is to a large degree defined by a purposeful neural process [for example, as a possible exploratory mechanism [START_REF] Roebuck | Neuromuscular responses to disturbance of balance in children with prenatal exposure to alcohol[END_REF]], this process may be suppressed if the controller is impaired in postural control. For example, patients with Parkinson's disease can demonstrate decreased postural sway while postural instability is one of the cardinal features of this disorder [START_REF] Romero | Changes in postural control with aging and Parkinson' s disease[END_REF].

Both interpretations given above consider, in a different degree, that the observed drop in the ML sway is a reflection of a purposeful strategy of the controller that suppresses the natural sway process that may by itself lead to loss of balance under the additional constraints. In addition, the results of the present study corroborate the idea that the control of the body sway in AP and ML directions during unperturbed upright stance are achieved by two distinct neuromuscular mechanisms that can be influenced by the outcome of each other's action.

In summary, the findings of this study indicate that the inverted pendulum model is unable to completely explain the processes involved in the control of the quiet upright stance in the anteriorposterior and medial-lateral directions.

Alexandrov, A. V., Frolov, A. A., Horak, F. B., Carlson-Kuhta, P., & Park, S. (2005). Feedback equilibrium control during human standing. Biological Cybernetics, 93, 309-322. Alexandrov, A. V., Frolov, A. A., & Massion, J. (2001). Biomechanical analysis of movement strategies in human forward trunk bending. II. Experimental study. Biological Cybernetics, 84, 435-443. Aramaki, Y., Nozaki, D., Masani, K., Sato, T., Nakazawa, K., & Yano, H. (2001). Reciprocal angular acceleration of the ankle and hip joints during quiet standing in humans. Experimental Brain Research, 136, 463-473. Baratto, L., Morasso, P. G., Re, C., & Spada, G. (2002). A new look at posturographic analysis in the clinical context: Sway-density versus other parameterization techniques. Motor Control, 6, 246-270. Table 1: Mean (± S.E.) of the COP, rambling, and trembling RMS and mean speed for four immobilizations and two visual conditions.

Figure captions

Caption to the Tables

Table 2: Summary of the effects of joint immobilization on standing balance.

  Six trials were performed in each immobilization condition, three trials with eyes open and three with eyes closed. The order in which each visual condition was presented in each trial was randomized across participants: a trial with open (close) eyes was always followed by a trial with eyes closed (open).

  depicts the averaged normalized RMS and mean speed values of the COP (left-hand side), and the rambling (middle) and trembling (right-hand side) displacements in the AP direction. INSERT FIGURE 4 HERE Effect of immobilization on the RMS The AP RMS of the COP and the rambling and trembling displacements averaged across joint immobilization and visual conditions are depicted in the upper panels of the Fig. 4. MANOVA revealed main effects of joint immobilization [Wilks' Lambda = .477, F(9, 61) = 2.4 p < .05, ² = .21] and vision [Wilks' Lambda = .07, F(3, 7) = 31, p > .001, ² = .93], but no interaction between those factors [Wilks' Lambda = .634, F(9, 61) = 1.39, p > .05, ² = .14] on the RMS in AP direction. Univariate analyses revealed a main effect of joint immobilization on the RMS of the trembling displacement [F(3, 27) = 4.29, p > .05, ² = .32], and marginal main effects on the RMS of the COP [F(3, 27) = 2.55 p = .076, ² = .22] and rambling [F(3, 27) = 2.49, p = .082, ² = .22] displacements. Post hoc tests with appropriated Bonferroni corrections showed that the RMS of the trembling displacement was larger when the knee, hips, and trunk were constrained (CT condition) as compared to when no constraint was imposed to the joints (NC condition). Although a marginal effect of joint immobilization on the RMS of the COP and rambling displacements was identified, no differences among joint immobilization conditions were observed. With respect to the main effect of vision, univariate analysis indicated that the RMS of the COP [F(1, 9) = 56, p < .001, ² = .86], and the rambling [F(1, 9) = 27, p < .005, ² = .75] and trembling [F(1, 9) = 48 p < .001, ² = .84] displacements were larger when participants kept their eyes closed (CE) as compared to when they kept their eyes open (OE).Effect of immobilization on the mean speedThe AP mean speed of the COP and the rambling and trembling displacements averaged across joint immobilization and visual conditions are depicted in the lower panels of Fig.4. For these variables, MANOVA revealed a marginal main effect of joint immobilization [Wilks' Lambda = .572, F(9, 61) = 1.75, p = .098, ² = .17] and a main effect of vision [Wilks' Lambda = .09, F(3, 7) = 24, p > .001, ² =.91], but no interaction between those factors [Wilks' Lambda = .72, F(9, 61) = .99, p > .05, ² = .10]. Univariate analyses revealed main effects of joint immobilization on mean speed of the COP [F(3, 27) = 3.77, p = .022, ² = .29] and trembling [F(3, 27) = 4.78, p = .008, ² = .35] displacements, but not on the mean speed of the rambling [F(3, 27) = .88, p = .462, ² = .21] displacement. For both, COP and trembling, post hoc tests revealed that the mean speed of the COP and the trembling displacements was larger in CT than in NC condition. Regarding the main effect of vision, univariate analysis indicated that the main speed of the COP [F(1, 9) = 70, p < .001, ² = .89], and the rambling [F(1, 9) = 70, p < .001, ² = .89], and trembling [F(1, 9) = 56, p < .001, ² = .86] displacements was larger in the CE than in OE condition.

  Fig. 5 depicts the averaged normalized RMS and mean speed values of the COP (left-hand side) and the rambling (middle) and trembling (right-hand side) displacements in the ML direction. INSERT FIGURE 5 HERE Effect of immobilization on the RMS The ML RMS of the COP and the rambling and trembling displacements averaged across immobilization and visual conditions are depicted in the upper panels of Fig. 5. MANOVA revealed significant main effects of joint immobilization [Wilks' Lambda = .346, F(9, 61) = 3.7, p < .005, ² = .30] and vision [Wilks' Lambda = .214, F(3, 7) = 8.55, p < .05, ² = .79] on RMS in ML direction. No interaction between joint immobilization and vision [Wilks' Lambda = .713, F(9,61) = 1, p > .05, ² = .41] was observed. Univariate analysis revealed a significant effect of joint immobilization on the RMS of the COP [F(3, 27) = 9.24, p < .001, ² = .51] and the rambling [F(3, 27) = 7.02, p < .005, ² = .43] and trembling [F(3, 27) = 10.6, p < .001, ² = .54] displacements. Post hoc tests revealed that the RMS of COP and the rambling and trembling displacements were smaller in the CH and CT than in the NC condition. Regarding vision, results revealed that RMS of the COP [F(1, 9) = 23.73, p < .005, ² = .73] and the rambling [F(1,9) = 17.2, p < .005, ² = .66] and trembling [F(1, 9) = 22.4, p < .005, ² = .71] displacements were larger in the CE than in the OE condition. Effect of immobilization on the mean speed The ML mean speed of the COP and the rambling and trembling displacements averaged across immobilization and visual conditions are depicted in the lower panels of the Fig. 5. MANOVA revealed significant effects of joint immobilization [Wilks' Lambda = .35, F(9, 61) = 3.66, p < .005, ² = .30], and vision [Wilks' Lambda = .271, F(3, 7) = 6.27, p < .05, ² = .73], and no interaction between those factors [Wilks' Lambda = .591, F(9, 61) = 1.63, p > .05, ² = .16] on the mean speed. Univariate analysis revealed significant effects of joint immobilization on the mean speed of the COP [F(3, 27) = 10.4, p < .005, ² = .51] and the rambling [F(3, 27) = 14.6, p < .001, ² = .43] and trembling [F(3, 27) = 7.1, p < .05, ² = .54] displacements. Post hoc tests revealed that the mean speed of the COP and the rambling and trembling displacements were smaller in the CH and CT than in the NC condition. Regarding vision, results revealed that the mean speed of the COP [F(1, 9) = 24, p < .005, ² = .73] and the rambling [F(1, 9) = 16, p < .005, ² = .64] and trembling [F(1, 9) = 11, p < .01, ² = .55] displacements were larger in the CE than in the OE condition.

Fig. 1 :

 1 Fig. 1: Schematic representation of the different immobilization conditions. No constraint (NC), constraint of the knee joints (CK), constraint of the knee and hip joints (CH), and constraint of the knees, hips, and trunk (CT).

Fig. 2 :

 2 Fig. 2: Sequence of operations used for COP trajectory decomposition. In A, the horizontal force and the instances when F hor = 0 (•) are presented. The COP displacement, the COP position at the IEP instants (•), rambling trajectory (estimated by the cubic spline interpolation of the COP position at the IEP instants), and trembling trajectory (the difference between COP and interpolated rambling trajectories) are presented in B. Participant SS, female, age 29 yrs; Ht 1.61 m; Wt 52 kg.

Fig. 3 :

 3 Fig. 3: A representative example of the results of the COP trajectory during all joints immobilization and visual conditions. Bars and numbers represent the standard deviation in the anterior-posterior (COP AP ) and medial-lateral (COP ML ) directions. Participant SS, female, age 29 yrs; Ht 1.61 m; Wt 52 kg.

Fig. 4 :

 4 Fig. 4: The RMS and speed values of the COP, rambling and trembling in the anterior-posterior direction averaged across participants. Values were normalized by the values of the NCOE condition (NCOE = 0%). Error bars indicate the standard error of the means. * represents differences with p < .05.

Fig. 5 :

 5 Fig. 5: The RMS and speed values of the COP, rambling and trembling in the medial-lateral direction averaged across participants. Values were normalized by the values of the NCOE condition (NCOE = 0%). Error bars indicate the standard error of the means. * represents differences with p < .05.

Fig.1

Fig.2

Fig.3

Fig.4

Fig.5

Table.1

Table.2