Corresponding author), U. Ivan Giorgio
email: ivan.giorgio@uniroma1.it

Ugo Andreaus

Daria Scerrato

Piero Braidotti

Modeling of a non-local stimulus for bone remodeling process under cyclic load: Application to a dental implant using a bioresorbable porous material

Keywords: Dental implant, cycle load, resorbable biomaterial, dissipation, visco-poro-elastic material, bone remodeling

may come

Modeling of a non-local stimulus for bone remodeling

process under cyclic load: Application to a dental implant using a bioresorbable porous material

Introduction

Dental implants are artificial systems usually consisting of an endosteal component, which is completely implanted in the mandible or jawbone and of an abutment which connects the endosseous component with the oral cavity, in order to replace one or more missing teeth [START_REF] Soncini | Computational approach for the mechanical reliability of a dental implant[END_REF]. The most recent experimental researches have been aimed at the study of the response of bone tissue following the insertion of an implant. In particular, researchers tried to make functional and structural interfaces between implant and bone and, then, to improve the process of osseointegration through surface treatments and design optimization of the system [START_REF] Andreaus | Coupling image processing and stress analysis for damage identification in a human premolar tooth[END_REF][START_REF] Ancillao | Finite element analysis of the stress state produced by an orthodontic skeletal anchorage system based on miniscrews[END_REF][START_REF] Marangos | Physico-mechanical properties determination using microscale homotopic measurements: application to sound and caries-a↵ected primary tooth dentin[END_REF]. The reaction of the bone tissue, in fact, depends on chemical properties of the materials used and on the ability of the implant to induce a mechanical state of stress in the bone tissue that can promote osteogenesis processes. Research, today, follows two main routes: i) clinical trial designed to evaluate the incidence of known failure mechanisms of the implant, and ii) preclinical study that allows to conduct a preliminary analysis of the system. This work follows the latter route: in fact, the use of numerical simulation based on the finite element method appears advantageous in terms of time and costs compared to carrying out experiments in vivo or in vitro on dental implants. In addition, it seems useful to guide the clinical trial with a gallery of examples obtained by numerical studies. In relation to issues of possible failure in the medium to long term, this paper is interested in the study of the remodeling of bone in the neighborhood of the metal stem as its resorption may cause disconnection of the fixture from the host tissue and therefore implant mobilization. The issues are addressed with the use of computational tools and are intended to define design criteria to increase the average life of an implant. In particular, in this paper, we numerically study an implant model characterized by axial symmetry and composed of fixture, trabecular bone tissue and porous biomaterial. It allows to determine the evolution of the states of deformation, stress, energy and dissipation that are originated in the compound. The geometry of the device that inspired this research is the one of the IntraMobil Zylinder (IMZ) and the Bio-Vent R implant [START_REF] Niznick | Product profile-achieving osseointegration in soft bone: The search for improved results[END_REF][START_REF] Babbush | Intramobile cylinder (IMZ) two-stage osteointegrated implant system with the intramobile element (IME): part I. Its rationale and procedure for use[END_REF][START_REF] Byrne | Fundamentals of Implant Dentistry[END_REF]. They are particularly suitable for soft and atrophic bone tissue, and are both characterized by coated cylinder implant designed for the mandibular jaw; the IMZ implant is coated with titanium plasma spray, the Bio-Vent implant is coated with hydroxyapatite.

Modeling

A dental implant consists of a titanium cylinder or screw, between 8 and 16 mm long, inserted in the jaw into a predisposed bony socket and acts as a replacement root for supporting a false tooth. Specifically, it is made of a titanium alloy, a material that is known to be compatible with biological tissues and is able to bond with adjacent bone during the healing process. Indeed at this stage, the titanium surface of the implant fuses with the surrounding bone, in a process known as osseointegration. Only after this stage, which can take about 3-6 months, the implant is stable enough to support a new tooth. The titanium alloy presents also advantageous mechanical properties, namely, strength, endurance and good fatigue strength.

In addition, this alloy can be fabricated with high precision to provide a wide variety of implants taking into account dimensions and state of health of bones and gums of a broad range of patients to ensure an optimal implant stability. Indeed, the problem of buckling due to the architecture of trabecular bone and the interaction with the graft under critical applied loads can occur. Possible strong problem singularities and instabilities can be regularized using methods described in the literature [START_REF] Di Egidio | Linear and non-linear interactions between static and dynamic bifurcations of damped planar beams[END_REF][START_REF] Luongo | On the e↵ect of twist angle on nonlinear galloping of suspended cables[END_REF][START_REF] Ruta | A beam model for the flexural-torsional buckling of thin-walled members with some applications[END_REF][START_REF] Rizzi | On the postbuckling analysis of thin-walled frames[END_REF][START_REF] Rizzi | Initial postbuckling behavior of thin-walled frames under mode interaction[END_REF].

Often there is a need to increase the quantity of bone at the jawbone site with a graft to allow an optimal connection of the implant. The reason of this might be, for instance a lack of bone at that site due to a long-term missing tooth, a periodontal disease and aging or injury. In these cases, it is possible to use a synthetic bone substitute (resorbable alloplast) that acts as a sort of sca↵old for the formation of new bone from the surrounding region and may be replaced, in part or completely, by natural bone.

In this paper, we consider a simple geometry to describe the interaction between a titanium insert, a bio-resorbable graft and the bone. In detail, we take into account a cylindrical region inside the jawbone which hosts the graft, shaped according to a half-ellipsoid, and the titanium insert (see Fig. 1). Cylindrical coordinates are the natural choice for this system exhibiting axial symmetry. Therefore, we introduce the usual coordinate for this reference frame characterized by the unit basis vectors (e r , e ✓ , e z ). Following a Lagrangian description we identify the constituent particles of the system by their position in a geometric configuration taken as a reference and denoted with C ⇤ , that is to say by the variable vector X whose coordinates in the adopted reference frame are (R, ⇥, Z). Then, we express the value of any physical value in the current configuration, i.e., C t , depending on the particle to which it is related and the current time, that is to say depending on variables X and t. In particular, a reference configuration is chosen such that the system is free-stress. However, since the system under study is characterized by a mass variable in time, the material properties of the reference configuration are not fixed, but evolve as the mass is varying. In addition, it is worth noting that the time scale of biological phenomena, namely, the evolution of the mass due to the remodeling process, is very slow compared to the scale characterizing the mechanical behavior and, therefore, the biological evolution of the system can be assumed not to entail the occurrence of any inertial e↵ect due to mass variation.

By assuming for axisymmetric reasons that all the considered fields are invariant under rotations, the position vector x = (r, 0, z) of the particle located in C ⇤ at position X is given by

x = re r + ze z = (R, Z, t) = Re r + Ze z + u(R, Z, t) (1)
where the displacement u of the particle located in C ⇤ at position X is obtained by

u(R, Z, t) = u r (R, Z, t)e r + u z (R, Z, t)e z (2)
We assume the transformation of the medium between the configurations C ⇤ and C t to be infinitesimal, namely, kru(X, t)k ⌧ 1, and therefore the Green-Lagrange strain tensor is reduced to the linearized strain tensor E that in cylindrical coordinates is expressed by the non-null components as

E rr = @u r @R , E zz = @u z @Z , E rz = 1 2 ✓ @u r @Z + @u z @R ◆ , E ✓✓ = u r R (3)
To properly describe the behavior of the porous materials involved in this system, namely, living bone and artificial bio-resorbable material, we consider a solid mixture composed of two phases (see e.g. [START_REF] Grillo | Evolution of a fibre-reinforced growing mixture[END_REF][START_REF] Placidi | Variational formulation of pre-stressed solid-fluid mixture theory, with an application to wave phenomena[END_REF][START_REF] Placidi | Thermodynamics of polycrystalline materials treated by the theory of mixtures with continuous diversity[END_REF]), endowed with a microstructure that takes the presence of connected pores into account. As a result, in the framework of micro-morphic continua (see e.g. [START_REF] Goda | A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization[END_REF][START_REF] Altenbach | Analysis of the viscoelastic behavior of plates made of functionally graded materials[END_REF][START_REF] Altenbach | Acceleration waves and ellipticity in thermoelastic micropolar media[END_REF][START_REF] Aminpour | A one-dimensional continuum with microstructure for single-wall carbon nanotubes bifurcation analysis[END_REF]), we introduce another kinematic descriptor able to characterize the mechanical response of pores inside the solid matrices considered (see [START_REF] Cowin | Linear elastic materials with voids[END_REF][START_REF] Biot | Mechanics of deformation and acoustic propagation in porous media[END_REF][START_REF] Biot | Generalized theory of acoustic propagation in porous dissipative media[END_REF][START_REF] Cowin | Bone poroelasticity[END_REF] and for further developments [START_REF] Dell'isola | A micro-structured continuum modelling compacting fluidsaturated grounds: The e↵ects of pore-size scale parameter[END_REF][START_REF] Dell'isola | Saint-Venant's problem for porous linear elastic materials[END_REF][START_REF] Madeo | A continuum model for deformable, second gradient porous media partially saturated with compressible fluids[END_REF][START_REF] Madeo | A variational deduction of second gradient poroelasticity II: An application to the consolidation problem[END_REF][START_REF] Tomic | Poroelastic materials reinforced by statistically oriented fibresnumerical implementation and application to articular cartilage[END_REF][START_REF] Federico | Possible approaches in modelling rearrangement in a microstructured material[END_REF][START_REF] Grillo | A poroplastic model of structural reorganisation in porous media of biomechanical interest[END_REF][START_REF] Dell'isola | Dynamics of solids with microperiodic nonconnected fluid inclusions[END_REF]): the change of the Lagrangian porosity ⇣, i.e. the change of the e↵ective volume of the fluid content per unit volume of the body, or mathematically

⇣(R, Z, t) = ( (R, Z, t), t) ⇤ (R, Z, t) (4)
where ⇤1 stands for the porosity in the reference configuration and is the Lagrangian porosity. The last can be expressed in terms of volume fraction of bone and artificial material, & b and & m respectively, as

= 1 (& b + & m ) = 1 ✓ ⇢ b ⇢b + ⇢ m ⇢m ◆ (5)
where ⇢ b , ⇢ m are the apparent mass densities of the two phases of the mixture, and ⇢b , ⇢m are the true mass densities (with self-explanatory indexes). This kind of approach indeed aims to incorporate microscale e↵ects, even if complex, into a continuum description (see e.g. [START_REF] Grillo | Evolution of a fibre-reinforced growing mixture[END_REF][START_REF] Gangho↵er | A contribution to the mechanics and thermodynamics of surface growth. Application to bone external remodeling[END_REF][START_REF] Gangho↵er | Mechanical modeling of growth considering domain variation-part II: Volumetric and surface growth involving eshelby tensors[END_REF][START_REF] Placidi | On the role of grain growth, recrystallization and polygonization in a continuum theory for anisotropic ice sheets[END_REF][START_REF] Placidi | An anisotropic flow law for incompressible polycrystalline materials[END_REF][START_REF] Dell'isola | What are the dominant thermomechanical processes in the basal sediment layer of large ice sheets?[END_REF] for applications to growth).

At the initial stage, the porous system is made of two distinct regions in which each of the phases, the artificial bio-resorbable material and the living bone, are separately present (see Fig. 1); it is only during the bone remodeling, in the process called osteoconduction, that new bone grows from the surrounding bone in the material region and therefore this zone exhibits both phases of the mixture considered. It is noted that the typical tools of the study of the phase transition can also be used in this context to model the rise of a new region in which both constituents coexist (see e.g. [START_REF] Dell'isola | A phenomenological approach to phase transition in classical field theory[END_REF][START_REF] Dell'isola | On phase transition layers in certain micro-damaged two-phase solids[END_REF][START_REF] Yeremeyev | The stability of the equilibrium of two-phase elastic solids[END_REF][START_REF] Eremeyev | Phase transitions in thermoelastic and thermoviscoelastic shells[END_REF]).

Following the same approach as proposed in [START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF], we assume the hypothesis of isotropic and inhomogeneous material; moreover, denoting respectively by ⇢ ⇤ b and ⇢ ⇤ m the mass densities of bone and of bio-resorbable material in the reference configuration, the stored energy density E (E, ⇣; ⇢ ⇤ b , ⇢ ⇤ m ) associated with strain tensor E and fluid volume distortion from the reference configuration ⇣ can be represented as

E = E s (E; ⇢ ⇤ b , ⇢ ⇤ m ) + 1 2 K 1 (⇢ ⇤ b , ⇢ ⇤ m )⇣ 2 + 1 2 K 2 kr⇣k 2 K 3 (⇢ ⇤ b , ⇢ ⇤ m ) ⇣ trE (6)
where

E s = 1 2 T el ij E ij
is the bulk strain-dependent energy density, and T el ij is the second Piola stress tensor given by the constitutive relation:

T el ij = (⇢ ⇤ b , ⇢ ⇤ m ) trE ij + 2µ(⇢ ⇤ b , ⇢ ⇤ m ) E ij (7)
1

From this point on, the superscript ⇤ denotes all quantities in the reference configuration.

being and µ the Lamé parameters and where we associate the values taken by the subscripts i and j ranging over the set {1, 2, 3}, in the given order, to the directions of the reference frame r, ✓, z in order to use the summation convention. Eq. ( 7) can be written in terms of Young's modulus, Y , and Poisson's ratio, ⌫, by recalling the following relationships:

= Y (⇢ ⇤ b , ⇢ ⇤ m )⌫ (1 2⌫) (1 + ⌫) µ = Y (⇢ ⇤ b , ⇢ ⇤ m ) 2 (1 + ⌫) (8)
As a first approximation, Poisson's ratio is assumed to be constant, while we hypothesize that Young's modulus changes with the bone remodeling according to the power-law:

Y = Y Max b (& b ⇤ ) b + Y Max m (& m ⇤ ) m (9)
where

Y Max b and Y Max m
are the maximal elastic moduli and the exponents b , m are constants. The coe cient K 1 which can be interpreted as a coe cient of compressibility related to the fluid phase inside pores, i.e. the bone marrow, is assumed to depend on the reference porosity as follows:

K 1 = ✓ ⇤ K f + (↵ B ⇤ )(1 ↵ B ) K dr ◆ 1 (10) where K f is the marrow modulus, K dr = Y /(3(1 2⌫
)) is the drained bulk modulus of the porous matrix and ↵ B is the Biot-Willis coe cient that satisfies the inequality ⇤ 6 ↵ B 6 1. Moreover, the trace of the small strain tensor

E trE = @u r @R + u r R + @u z @Z (11)
is introduced to describe the interaction between the solid and the microstructure. The gradient of the change of porosity, because of the axial symmetry, is given by

r⇣(R, Z, t) = @⇣ @R e r + @⇣ @Z e z ( 12
)
and is used as an application of a second gradient theory -Cahn-Hilliard fluids-to take capillary interaction phenomena among neighbouring pores into account [START_REF] Seppecher | Second-gradient theory: application to Cahn-Hilliard fluids[END_REF]. The material parameter K 2 that appears in Eq. ( 6) can be thus interpreted as a sti↵ness related to the influence of capillarity governed by the r⇣. Higher order theories of deformation can also be usefully employed to model damage evolution of bone (see e.g. [START_REF] Yang | Higher-order stress-strain theory for damage modeling implemented in an element-free galerkin formulation[END_REF][START_REF] Yang | Higher-order continuum theory applied to fracture simulation of nanoscale intergranular glassy film[END_REF][START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF][START_REF] Placidi | A variational approach for a nonlinear one-dimensional damage-elasto-plastic secondgradient continuum model[END_REF][START_REF] Rinaldi | A microscale second gradient approximation of the damage parameter of quasi-brittle heterogeneous lattices[END_REF][START_REF] Andreaus | Damage modelling and seismic response of simple degrading systems[END_REF][START_REF] Carcaterra | Macroscopic description of microscopically strongly inhomogenous systems: A mathematical basis for the synthesis of higher gradients metamaterials[END_REF]) which may be relevant in the remodeling process [START_REF] Lee | Bone adaptation to load: microdamage as a stimulus for bone remodelling[END_REF].

The material parameter K 3 introduced in Eq. ( 6) can be physically interpreted as the coupling between microstructure and solid bulk [START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF] and is evaluated according to the relation:

K 3 = p ĝ( ⇤ ) K 1 ( 13
)
where the monotonic function ĝ( ⇤ ) takes non-negative value less then one (see Fig. 2) and is designed to emphasize the e↵ect of pores when the level of porosity is high and to understate this e↵ect when porosity decreases according to: where A k 3 and s k 3 are shape coe cients that can be chosen to characterize the coupling law adopted. Specifically, the coe cient A k 3 modulates the amplitude of the coupling law Eq. ( 14) while s k 3 governs the slope of the above mentioned law. Recalling that the material parameter K 3 governs the coupling between the variables of the solid bulk and those of the microstructure, the coupling law Eq. ( 14) has been designed in order to i) penalize the coupling when the porosity tends to zero, and then the global behavior tends to the one of the solid bulk, and ii) exalt the coupling when the porosity tends to one, and then the role of the microstructure becomes increasingly more important.

ĝ( ⇤ ) = A k 3 ⇡ ⇢ atan  s k 3 ✓ ⇤ 1 2 ◆ + atan ⇣ s k 3 2 ⌘ 0 < A k 3 6 1 ( 14
)
In order to guarantee the positive definiteness of the stored energy density (6), the following inequalities which involve the material parameters above defined should be satisfied as stated in [START_REF] Cowin | Linear elastic materials with voids[END_REF]:

µ > 0,  = + 2/3µ > 0, K 1 > 0, K 2 > 0, K 1 > 4K 2 3 ( 15
)
Finally, interface conditions for the variable ⇣ can be formulated by adding to the energy density of Eq. ( 6) a boundary extra term:

E int = 1 2 K 4 ⇣ + ⇣ 2 = 1 2 K 4 [[⇣]] 2 with K 4 > 0 ( 16
)
which accounts for the contribution of the interface, @B int , between bone and material regions. Therefore, the material parameter K 4 plays the role of a sti↵ness that rules the interaction between the two materials in terms of the jump of ⇣. For the sake of simplicity, the coe cients K 2 and K 4 are assumed to be constant. Some features suggests the presence of viscous damping in bone [START_REF] Sedlin | A rheologic model for cortical bone: a study of the physical properties of human femoral samples[END_REF]: total deformation under a constant stress was not attained instantaneously, i.e, bone has a retardation time; an hysteresis loop arises from loading-unloading cyclic tests; and so forth. Therefore, according to the Kelvin-Voigt formulation, viscous damping can be added to the material model of the porous mixture by means of an extra stress T v ij related to the rate of elastic strain of the material

T v ij ( Ė) = 2µ v ✓ Ėij 1 3 Ėkk ij ◆ +  v Ėkk ij (17)
where  v and µ v are the bulk and shear viscosity coe cients, respectively. As a result, the dissipated power is given by 2D

s = T v ij Ėij (18)
It is worth observing that the validity of this simple model, satisfactory for low physiological rate evolutions [START_REF] Bargren | Mechanical properties of hydrated cortical bone[END_REF], becomes increasingly inadequate beyond the range of quasi-static loading which is herein considered (for more complex models see e.g. [START_REF] Carcaterra | Dissipation in a finite-size bath[END_REF][START_REF] Carcaterra | Fractional dissipation generated by hidden wave-fields[END_REF]).

Since the problem under study is time-dependent, but the evolution in time is very slow, a quasi-static formulation is assumed and then inertial e↵ects are neglected. The numerical code used to calculate the configuration sequences of the considered system includes a formulation of the poro-viscoelastic problem linking external mechanical applied loads to the corresponding equilibrium deformation states. However, it is also required the knowledge of the evolution of the mass densities, from which the sti↵ness of the mixture, given by Eq. ( 9), can be evaluated at each step. The evolution problem involves only the equations for time variations of apparent mass densities of the living bone tissue and the bio-resorbable material. On the other hand, since time variations of mixture sti↵ness and variations of loads occur in a time scale which is much larger than the transient characteristic time leading to macroscopic mechanical equilibrium, the two problems, the mechanical and the evolutionary one, can be separately solved at each time step in sequence with an iterative scheme as time increases.

Poro-viscoelastic formulation

The mechanical governing equations of the graft/bone system can be deduced from the generalized principle of virtual work. Indeed, the total virtual work, including dissipative actions (see e.g. [START_REF] Lekszycki | Application of variational methods in analysis and synthesis of viscoelastic continuous systems[END_REF]), for any arbitrary virtual displacement u and for a virtual change of porosity ⇣ equal to zero is given by

Z B h T el ij E ij + K 1 ⇣ ⇣ + K 2 ⇣ ,i ⇣ ,i K 3 (E ii ⇣ + ⇣ E ii ) i dV Z B T v ij E ij dV Z @B int K 4 [[⇣]] [[⇣]]dS + Z @⌧ B ⌧ i u i dS + Z @B ⌅ ⇣dS = 0 (19)
To solve our problem, we consider a mixed boundary value problem in which we impose the displacement components

u r (R e , Z, t) = 0 8 Z and t (20) u z (R, 0, t) = 0 8 R and t ( 21
)
being R e the external radius of the cylinder under study. The titanium stem (here modeled as a linear elastic solid) is loaded by an external pressure on its top, i.e. @ ⌧ B T (see Fig. 1); in turn, the titanium insert applies a surface traction ⌧ i on the boundary @ ⌧ B between the titanium stem itself and the mixture. Possibly a micro-structural action ⌅, related to the dilatant behavior of the porous mixture, can be applied on the boundary @B. For further information about contact interactions in more general cases see for example [START_REF] Dell'isola | Cauchy tetrahedron argument applied to higher contact interactions[END_REF][START_REF] Dell'isola | The postulations á la D'Alembert and á la Cauchy for higher gradient continuum theories are equivalent: a review of existing results[END_REF]. Since inertia e↵ects are neglected but the mechanical governing equation is a di↵erential equation of first order in time for the presence of the dissipative actions, it is worth recalling that herein the initial conditions are reduced to

u(R, Z, 0) = 0 8 R, Z (22)

Evolution rules for the mass densities

The evolution equations for apparent mass densities are assumed to be first order ordinary di↵erential equations with respect to time. We also assume that the spatial gradient of mass densities does not a↵ect the remodeling equations, and therefore we simply assume that the process of remodeling is non-local only as a function of the nature of the stimulus (herein it is obtained by an integral operator), and it does not depends on what is occurring outside the representative elementary volume by other non-local phenomena of pure mechanical nature. As a result, we assume that the total mass density rates depend on the stimulus and the porosity field as follows:

( ⇢⇤ b = A b (S) H ( ) with 0 < ⇢ ⇤ b 6 ⇢b ⇢⇤ m = A m (S) H ( ) with 0 < ⇢ ⇤ m 6 ⇢ 0 m ( 23
)
where A b and A m are piece-wise linear functions with di↵erent slopes for negative (r b and r m ) and positive values of stimulus (s b ); moreover for the bio-resorbable material and positive values of the stimulus the function A m vanishes, since synthesis of bio-resorbable material cannot occur. Indeed, the maximum value possible for ⇢ ⇤ m is its initial one, i.e. ⇢ 0 m . Specifically, we have

A b (S) = ⇢ s b S for S 0 r b S for S < 0 (24) A m (S) = ⇢ 0 for S 0 r m S for S < 0 (25)
The function H (see Fig. 3) is evaluated as a normalized Specific Surface, i.e., the internal surface area per unit volume of whole bone [START_REF] Giorgio | A visco-poroelastic model of functional adaptation in bones reconstructed with bio-resorbable materials[END_REF]. Indeed, this function is related to the surface available for resorption or synthesis 'inside' the considered material particle which, in its turn, depends on the 'e↵ective' porosity field [START_REF] Lekszycki | A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials[END_REF].

To take the experimental evidence into account, we consider a so called 'lazy zone' for stimulus. Particularly, when the signal falls in a given range then the e↵ect of the actor cells is not macroscopically appreciable. For this purpose, we introduce two activation thresholds, P s ref and P r ref ('s' stands for synthesis while 'r' means resorption), that define the width of this 'lazy zone' associated with a balanced biological equilibrium state and, therefore, the stimulus becomes: where P (X, t) is

S(X, t) = 8 > < > : P (X,
P (X, t) = R B [a E s (X 0 , t) + b D s (X 0 , t)] $ [⇢ b (X 0 , t)] e kX X 0 k 2 2D 2 dX 0 R B e kX X 0 k 2 2D 2 dX 0 (27)
We assume that the signal stemming from sensor cells in a given material particle of the bone/bio-material mixture is instantaneously transmitted, i.e., the transmission time scale is negligible compared to the characteristic time of the remodeling phenomena, and that its intensity decays exponentially with the distance from it [START_REF] Lekszycki | A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials[END_REF][START_REF] Mullender | Proposal for the regulatory mechanism of wol↵'s law[END_REF][START_REF] Madeo | A continuum model for the bio-mechanical interactions between living tissue and bio-resorbable graft after bone reconstructive surgery[END_REF]. Besides, to emulate the non-local behavior due to the sensor cell network, we express the stimulus as a spatially averaged signal over a spherical region of radius D which represents the range of action of sensor cells [START_REF] Kumar | Dissipation energy as a stimulus for cortical bone adaptation[END_REF]. The function $ is assumed to be:

$ (⇢ ⇤ b ) = ⌘ tanh(⇠ ⇢ ⇤ b ) 0 < ⌘ 6 1 (28)
and it can be interpreted, in the framework of the feedback control theory [START_REF] Turner | Homeostatic control of bone structure: an application of feedback theory[END_REF], as a gain for the actual mechanical stimulus, i.e., 'the actuating signal', that takes the real activities of the osteocytes into account [START_REF] Giorgio | A visco-poroelastic model of functional adaptation in bones reconstructed with bio-resorbable materials[END_REF]. Several stimulus models have been proposed in literature to match experimental observations. All these models are characterized by a particular description of the stimulus; in some of them it is based on the strain energy density [START_REF] Lekszycki | A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials[END_REF][START_REF] Mullender | Proposal for the regulatory mechanism of wol↵'s law[END_REF][START_REF] Weinans | The behavior of adaptive bone-remodeling simulation models[END_REF], on a strain-damage energy density [START_REF] Hambli | Connecting mechanics and bone cell activities in the bone remodeling process: an integrated finite element modeling[END_REF], on the adaptation of the elastic modulus of the bone coupled directly to deviations of the strain tensor [START_REF] Cowin | Bone remodeling I: theory of adaptive elasticity[END_REF], while others propose a creepfatigue cumulative damage model [START_REF] Carter | A cumulative damage model for bone fracture[END_REF][START_REF] Caler | Bone creep-fatigue damage accumulation[END_REF], or a global optimal criterion for remodeling [START_REF] Lekszycki | Modelling of bone adaptation based on an optimal response hypothesis[END_REF][START_REF] Jang | Application of design space optimization to bone remodeling simulation of trabecular architecture in human proximal femur for higher computational e ciency[END_REF][START_REF] Andreaus | Optimal-tuning pid control of adaptive materials for structural e ciency[END_REF][START_REF] Andreaus | An optimal control procedure for bone adaptation under mechanical stimulus[END_REF][START_REF] Andreaus | Modeling of trabecular architecture as result of an optimal control procedure[END_REF][START_REF] Andreaus | Optimal bone density distributions: Numerical analysis of the osteocyte spatial influence in bone remodeling[END_REF]. In this paper, we propose for bone adaptation the use of a linear combination of the strain-dependent energy, E s , and the dissipation power, D s , in order to predict two important features of the remodeling phenomenon: the bone capability to control its mass and sti↵ness directly related to static mechanical demands [START_REF] Lekszycki | A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials[END_REF][START_REF] Weinans | The behavior of adaptive bone-remodeling simulation models[END_REF] for the function of load-carrying capacity, and the e↵ects of the frequency of applied loads [START_REF] Lanyon | Static vs dynamic loads as an influence on bone remodelling[END_REF][START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF] for the function of damping and shock absorption.

As a note, herein, we consider internal remodeling as being the resorption or deposition of bone material, accompanied by the removal and densification of the architecture of bone, but without any change in the overall shape [START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF][START_REF] Lekszycki | A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials[END_REF][START_REF] Andreaus | A 2-D continuum model of a mixture of bone tissue and bio-resorbable material for simulating mass density redistribution under load slowly variable in time[END_REF][START_REF] Giorgio | The influence of di↵erent loads on the remodeling process of a bone and bioresorbable material mixture with voids[END_REF]. An external or surface modeling is not studied. Such a kind of remodeling causes a change of the external shape of the overall bone structure, and occurs by the resorption or deposition of bone material on the surfaces (see e.g. [START_REF] Gangho↵er | A contribution to the mechanics and thermodynamics of surface growth. Application to bone external remodeling[END_REF][START_REF] Gangho↵er | Mechanical modeling of growth considering domain variation-part II: Volumetric and surface growth involving eshelby tensors[END_REF]).

The bio-material which we assume to use as a graft is absolutely not active in both sensing and actuating process. However, one could conceive a 'smart' prosthesis having much more e↵ective performances by using the concept developed in [START_REF] Giorgio | Piezo-electromechanical smart materials with distributed arrays of piezoelectric transducers: current and upcoming applications[END_REF][START_REF] Giorgio | Multimode vibration control using several piezoelectric transducers shunted with a multiterminal network[END_REF][START_REF] Porfiri | Modeling and design of passive electric networks interconnecting piezoelectric transducers for distributed vibration control[END_REF][START_REF] Maurini | Comparison of piezoelectronic networks acting as distributed vibration absorbers[END_REF][START_REF] Pagnini | The three-hinged arch as an example of piezomechanic passive controlled structure[END_REF] or some new metamaterials [START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF] which exhibit interesting features such as lightness and resistance to damage as proposed in [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF][START_REF] Dell'isola | Plane bias extension test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange multipliers and a perturbation solution[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF][START_REF] Madeo | Wave propagation in pantographic 2D lattices with internal discontinuities[END_REF].

Numerical simulations

In order to test the predictive capabilities of the proposed model, we performed numerical simulations, using the commercial software COMSOL Multiphysics R , on a three-dimensional cylindrical sample of height h b = 10 mm and radius R e = 2 h b /3. The sample is made of a piece of jawbone, a bio-material graft and a titanium insert as sketched in Fig. 1.

Initially, the mass density distributions of the bone and of the graft material are taken spatially uniform and segregated in two distinct zones, C b and C m respectively, to deal with a simple but representative case for a dental implant. Hence, initial conditions for bone mass density, ⇢ ⇤ b , and bio-material density,

⇢ ⇤ m , are ( ⇢ ⇤ b (X, 0) = ⇢ 0 b = ⇢b /2 8X 2 C b ⇢ ⇤ m (X, 0) = ⇢ 0 m = ⇢m /2 8X 2 C m ( 29
)
Herein, the problem under study is recast in a non-dimensional form by normalizing the stored elastic energy with respect to the maximal sti↵ness of bone Y Max b = 18 GPa and the dissipated power by Y Max b /t ref , being t ref a reference time assumed equal to 7.2576⇥10 6 s; the lengths are normalized with respect to the diameter of the sample, and the apparent mass densities with respect to the ⇢b = 1, 800 kg/m 3 , which is assumed equal to ⇢m . The non-dimensional quantities are denoted by a tilde, and can be summarized as:

Kf = K f /Y Max b , K2 = K 2 /(Y Max b 4R 2 e ), ˜ = /Y Max b , K4 = K 4 /(Y Max b 2R e ), μv = µ v /(t ref Y Max b ), v =  v /(t ref Y Max b ) (30)
A cyclic compression load case is considered, in which an external force is applied on the titanium insert along its longitudinal axis; in particular, the force ⌧ ext z (t) = F 0 F 1 sin(2⇡f t) with a pre-load F 0 and a cyclic contribution of magnitude F 1 = F 0 /2 and frequency f = {5, 10, 20} cycles per unit of time is . In absolute terms, the external force applied on the top of the fixture corresponds to a typical masticatory load of about 160 N on the abutment top. Young's modulus of titanium insert is assumed to be Y T = 108 GPa. The limit values that define the range of the 'lazy zone' are P r ref = 1.63 ⇥ 10 9 and P s ref = 1.84 ⇥ 10 9 . The values of the parameters used in the performed numerical simulations are summarized in Tabs. 1-3. Plots of some results relatively to the probe point P m ⌘ (R e /8, 10/13R e ) in the time range considered for simulation are showed since this point is located in the most stressed area, and it involves the evolution of both components of the mixture.

As a side note, recent applications of the finite element method (as discussed in [START_REF] Cazzani | Isogeometric analysis of plane-curved beams[END_REF][START_REF] Turco | A three-dimensional b-spline boundary element[END_REF][START_REF] Cazzani | Isogeometric analysis: a powerful numerical tool for the elastic analysis of historical masonry arches[END_REF][START_REF] Cazzani | Constitutive models for strongly curved beams in the frame of isogeometric analysis[END_REF][START_REF] Della Corte | Referential description of the evolution of a 2D swarm of robots interacting with the closer neighbors: Perspectives of continuum modeling via higher gradient continua[END_REF][START_REF] Federico | An energetic approach to the analysis of anisotropic hyperelastic materials[END_REF][START_REF] Solari | A numerical algorithm for the aerodynamic identification of structures[END_REF][START_REF] Greco | B-Spline interpolation of Kirchho↵-Love space rods[END_REF][START_REF] Greco | An implicit G1 multi patch B-spline interpolation for Kirchho↵-Love space rod[END_REF][START_REF] Greco | An isogeometric implicit G1 mixed finite element for Kirchho↵ space rods[END_REF][START_REF] Cuomo | A variational model based on isogeometric interpolation for the analysis of cracked bodies[END_REF]) could be usefully employed in this context to enhance the possibilities to improve the formulation and thus to ensure a better fit for experimental evidence.

Figure 4 shows the reference porosity ⇤ at the end of the process for di↵erent frequencies of loading. With reference to the region of the biomaterial, it is observed that for low frequencies, the process is slower and the biomaterial has more time to be reabsorbed. The area with high porosity (0.7 -blue zone) compared to initial value (0.5) is very wide at frequency of 5 cycles per unit of time and is gradually being reduced with increasing frequency; in this area the porosity of the mixture coincides with the one of the biomaterial; in the red area (which is the most stressed because it is below the implant) where the porosity decreases, there is growth of bone tissue and therefore the porosity is equal to the one of the mixture. For high frequencies, the growth of bone tissue invades the entire region of the biomaterial (red zone) and the minimum porosity (0.2) is reached.

Figure 5 shows the final mass densities of biomaterial represented in terms of variables R, Z at the end of the simulation for di↵erent frequencies. The figure indicates that in the region of the biomaterial (or at least in the more external discharged part) the reduction of the mass density of the material, in front of bone growth, decreases as the frequency increases (from 0.5 to 0.3 at 5 cycles per unit of time, 0.5 to 0.32 at 10 cycles per unit of time, and finally from 0.5 to 0.41 at 20 cycles per unit of time), because the bone tissue colonizes the material quicker as the frequency increases.

Figure 6 shows the change of the porosity ⇣ at the end of the process for di↵erent frequencies. The increase of frequency causes the mixture to exhibit a lower porous shrinkage leading the bone to grow faster and the pores to be filled quicker. Hence, the system will become more compact.

Figure 7 shows the time histories of the mass densities of bone (solid line) and material (dotted-dashed line) at the point P m probe for di↵erent frequencies. It is noted that while the mass density of the bone also grows significantly as the frequency increases, the mass density of the material remains almost constant and independent of frequency, as the bone is filling the voids of biomaterial.

Figure 8 shows the time histories of the change of porosity ⇣ at the probe point P m . As the frequency increases the mass density of mixture grows and the porosity decreases (see Fig. 4), therefore the mixture becomes sti↵er and hence the change of porosity ⇣ oscillates within an increasingly narrower band. The trend is clearly stabilizing towards the stationary state confirming the adequate choice of the analysis duration.

Figures 9 and10 show the time histories of the elastic energy E s and of the dissipation power D s at the probe point P m . It confirms the fact that, as the frequency increases, the mixture exhibits a less deformable and more dissipative behavior.

Figure 11 shows the energy dissipation loops for the three principal directions at the probe point P whose coordinates are (0.3R e , 1.154R e ). The diagrams were obtained by plotting the principal stresses in terms of the corresponding principal strains for di↵erent frequencies. The rows correspond to the principal directions and the columns to the loading frequencies. The area subtended by each cycle is proportional to the energy dissipated in the same cycle. The first cycle is marked with green color, the last cycle with red color. The shape of each cycle is approximately elliptical, due to the dissipation of a viscous nature. The problem is complex, depending on such factors that influence each other resulting interdependent as the stimulus, the dissipation, the evolution of the mass density of bone material which in turn a↵ects the sti↵ness of the mixture and thus the mechanical behavior. Starting from the third principal direction (third row), we note that both compressive stress and strain are negative. At low frequencies (5 cycles per unit of time), the deformation increases, the stress decreases and the final cycle is slightly larger than the initial one, that is, the dissipation grows. The growth process is slower; because less bone mass is produced and an erosion occurs in the area of the material (Fig. 4a) a higher porosity, a lower sti↵ness and a larger deformation are therefore attained. At high frequencies (20 cycles per unit of time), it is observed a di↵erent behavior compared to the one at low frequencies, by virtue of the faster evolution of bone mass density. In fact, deformation and dissipation significantly decrease, whereas a moderate stress reduction is exhibited. For intermediate frequencies, the two conflicting behaviors are approximately balanced and then final and initial cycles do not significantly di↵er in terms of both deformation and dissipation. Turning to the first principal direction (first row), characterized by deformation of elongation and tensile stress, it is observed a behavior similar to that already noted for the third principal direction. In fact, the factor which drives the described phenomenon is the rate of growth of the bone mass related to the frequency of the applied action. For low frequencies, deformation and dissipation increase while the stress decreases. For high frequencies, deformation, stress and dissipation decrease. For intermediate frequencies, the deformation and dissipation are approximately comparable. Finally, the second principal direction (second row) refers to the circumferential direction ✓. From a quantitative point of view, the deformations are of an order of magnitude smaller compared to the deformations in the other two principal directions (Poisson's e↵ect), but the overall behavior reflects the considerations already made with regard to the two other principal directions.

Figure 12 shows the influence of the width of the lazy zone on the reference porosity ⇤ of the mixture at the end of the process, at the frequency of 5 cycles per unit of time. The resorption threshold P r ref has been fixed and the synthesis threshold P s ref is increased from the value equal to P r ref .

It is observed that as the synthesis threshold increases, the red area characterized by low porosity (high mass density) visibly decreases in the most stressed area below the implant.

Conclusions

The subject of this work is a type of dental implant IMZ that is suitable to be inserted in soft and atrophic bone tissue of mandibular jaw. It was studied the interaction of the implant with the biomaterial and the functional adaptation of the bone tissue under a load slowly variable in time with di↵erent frequencies, that produces a mechanical non-local stimulus consisting in the linear combination of deformation energy and viscous dissipation. The main mechanical parameters that have been considered are: mass density, porosity, change of porosity, density of strain energy and viscous dissipation. Attention was focused on the evaluation of the influence of the load frequency and of the amplitude of the lazy zone. The increase of frequency within a suitable range favors the production of bone as already reported in many studies (see e.g. [START_REF] Lanyon | Static vs dynamic loads as an influence on bone remodelling[END_REF][START_REF] Turner | Three rules for bone adaptation to mechanical stimuli[END_REF]). However, if the stimulus is such to induce a too fast process of bone growth, then the conditions that should facilitate resorption of the biomaterial are not experienced. Therefore, the problem is to find a fair compromise between osteointegration and resorption of the biomaterial. The aim of this study was precisely to investigate the ways in which this compromise occurs and to quantify the limits of implementation. As regards the influence of the lazy zone, the performed numerical simulations have confirmed and quantified what might be expected intuitively, namely that the phenomenon of growth of the bone tissue is slowed down by an increase in the threshold of the stimulus. The numerical results presented in this study can be considered in view of reducing the risk of mobilization of the metallic stem with respect of the host bone tissue in the medium and long term; indeed, such a risk is due to the possibility of a non-perfect osteointegration at the interface between the implant cylinder and the jawbone.

Figure 1 :

 1 Figure 1: Initial arrangement of sample under study.

Figure 2 :

 2 Figure 2: Function ĝ( ⇤ ).

Table 2 :

 2 Values of the micro-parameters used in numerical simulations. non-dimensional value of the force is F0 = 0.00286 obtained normalizing by the sti↵ness Y Max b

  (a) frequency of 5 cycles per unit of time (b) frequency of 10 cycles per unit of time (c) frequency of 20 cycles per unit of time

Figure 4 :

 4 Figure 4: Reference porosity ⇤ at the end of the simulation.

c

  (a) frequency of 5 cycles per unit of time (b) frequency of 10 cycles per unit of time (c) frequency of 20 cycles per unit of time

Figure 5 :

 5 Figure 5: Normalized mass density of bio-material at the end of the simulation.

(a) frequency of 5

 5 cycles per unit of time (b) frequency of 10 cycles per unit of time (c) frequency of 20 cycles per unit of time

Figure 6 :

 6 Figure 6: Change of porosity ⇣ at the end of the simulation.

Figure 7 :

 7 Figure 7: Mass density histories of bone (solid line) and of bio-resorbable material (dotted-dashed line) at the probe point P m for di↵erent load frequencies (5, 10 and 20 cycles per unit of time).

Figure 8 :

 8 Figure 8: History of the change of the porosity ⇣ at the probe point P m for di↵erent load frequencies (5, 10 and 20 cycles per unit of time).

Figure 9 :

 9 Figure 9: History of the elastic energy at the probe point P m for di↵erent load frequencies (5, 10 and 20 cycles per unit of time).

Figure 10 :

 10 Figure 10: History of the dissipation at the probe point P m for di↵erent load frequencies (5, 10 and 20 cycles per unit of time).

  (a) freq. 5 cycles per unit of time (b) freq. 10 cycles per unit of time (c) freq. 20 cycles per unit of time (d) freq. 5 cycles per unit of time (e) freq. 10 cycles per unit of time (f) freq. 20 cycles per unit of time (g) freq. 5 cycles per unit of time (h) freq. 10 cycles per unit of time (i) freq. 20 cycles per unit of time

Figure 11 :

 11 Figure 11: Energy dissipation loops for the three principal directions ((a),(b),(c) -first principal direction; (d),(e),(f) -second principal direction; (g),(h),(i) -third principal direction) at the probe point P c .

9 Figure 12 :

 912 Figure 12: Influence of the lazy zone: reference porosity ⇤ at the end of the simulation and with a frequency load of 5 cycles per unit of time.

  t) P s

		1
		0.9
		0.8
		0.7
	function H	0.4 0.5 0.6
		0.3
		0.2
		0.1
		0	0	0.2	0.4	0.6	0.8	1
					porosity ?
					Figure 3: Function H.
					ref	for P (X, t) > P s ref
			0			for P r ref 6 P (X, t) 6 P s ref	(26)
			P (X, t) P r ref	for P (X, t) < P r ref
					8

Table 1 :

 1 Values of the macro-parameters used in numerical simulations.

Table 3 :

 3 Rates of resorption and synthesis.

	sb rb rm
	12 10 15

Acknowledgments

The research was funded by the Italian Ministry of University and Research, under the Scientific Research Program of Relevant National Interest: Year 2010-2011, Protocol 2010MBJK5B-005, Title "Dynamics, Stability and Control of Flexible Structures".


domain: phys spi

text:
published or not. The documents may come

MIMO radar transmit array fed by a 1 x M passive chaotic cavity

Ettien Lazare Kpré, Thomas Fromenteze, Cyril Decroze

MIMO radar transmit array fed by a 1 × M passive chaotic cavity Ettien L. Kpré 1 , Thomas Fromenteze 1 , Cyril Decroze 1 , 1 XLIM Research Institute -CNRS 123, Avenue Albert Thomas, 87060 Limoges Cedex, France e-mail : ettien.kpre@etu.unilim.fr

Abstract-It is well known that MIMO radar waveforms should be orthogonal to probe simultaneously the MIMO channel matrix H. This method is not easily implementable since it requires as much transmitter chains as antennas. Herein, a compressive technique is introduced to measure the channel information avoiding multiple waveforms generation. This technique is thus based on the generation of a single waveform to probe the channel by the mean of a passive chaotic component connected to the transmit array. Mathematical formulations are presented considering a device with uncorrelated transfer functions allowing the transmission of orthogonal waveforms. Simulation and experimental results are presented to highlight the advantages of such a technique. Index Terms-MIMO Radar, passive chaotic microwave component.

I. INTRODUCTION

MIMO radar became more popular and attracted much attention thanks to its ability to enhance radar performance [START_REF] Bliss | Multiple-input multiple-output (MIMO) radar and imaging: degrees of freedom and resolution[END_REF]. The essence of this concept is the possibility to probe simultaneously the channel with M orthogonal signals and record the backscattered signals with N receiving antennas as depicted in Fig. 1. The received signal can thus be reassigned to the sources by performing correlation in pairs. The channel matrix H M×N represents the interaction with targets measured at each frequency, required for the image reconstruction [START_REF] Jian | MIMO radar signal processing[END_REF]. Hence, each transmitting antenna should be addressed by a single waveform generated by a transmitter chain imposing the use of sophisticated and redundant RF systems. In this article, a single input MIMO radar is proposed to reduce the number of transmitters maintaining the same number of antennas. The principle is based on the use of a chaotic microwave component in which transfer functions are uncorrelated. The first approach introduced in [START_REF] Carsenat | UWB Antennas Beamforming Using Passive Time-Reversal Device[END_REF] was based a time reversal technique applied to a 1 × 4 small reverberation cell beamformer connected to UWB antennas. With a single input waveform, the authors demonstrated the control of antenna weightings without the need of active components or moving parts. This paper proposes the extension of this technique to MIMO radar transmit arrays. In this context, a single waveform is generated at the input of a chaotic component, and each transmit antenna can be addressed independently leading to the reduction of the number of transmitters compared to the conventional MIMO radar architectures. It is shown in this paper that the channel matrix can thus be retrieved from the chaotic component's transfer function matrix inversion. The remainder of this paper is organized as follows. Section II describes the conventional MIMO radar principle. Section III introduces the principle of unique input channel MIMO radar and the mathematical formulations ensuing. Simulation and experimental results are presented in section IV and V to present the advantages of the proposed method. Finally, section VI concludes this paper.

II. CONVENTIONAL MIMO RADAR

Let us consider a radar architecture as depicted in Fig. 1. Ideally, the radiating aperture array consists of M transmitting antennas radiating orthogonal waveforms and N receiving antennas. However, in practice, the generation of orthogonal waveforms can hardly be achieved by real-life hardwares. Numerous papers have thus been focused on MIMO radar pseudo-orthogonal waveform design [START_REF] Wang | MIMO SAR OFDM Chirp Waveform Diversity Design With Random Matrix Modulation[END_REF]- [START_REF] Babur | Nearly Orthogonal Waveforms for MIMO FMCW Radar[END_REF]. For example, the waveforms considered in this paper are frequency hopping modulated signals that can be computed considering Eq. ( 1):

s(t) = Q-1 ∑ q=0 u(t -qT c). exp( j2π( f 0 + c(q)∆ f )t) (1)
where:

u(t) = 1, if 0 < t < T c. 0, otherwise. ( 2
)
and c(q) is the q th element scrambled by a pseudo-noise Gold sequences [START_REF] Buehrer | Code Division Multiple Access (CDMA)[END_REF], and Q, the number of frequencies. T c corresponds to the pulse duration and is fixed by the total bandwidth Bw of the signal. The echo signals can be reassigned to the sources by performing correlation in pairs. The aim of the post processing is to estimate the channel matrix [H( f )] M×N which represents the interaction between the transmitting and receiving antennas:

H =    h 11 h 12 . . . h 1N h 21 h 22 . . . h 2N . . . . . . . . . . . . h M1 h M2 . . . h MN    (3)
where the channel between the antennas m and n is :

h mn ∝ exp( jk(x T m + x Rn )d sin θ 0 ). (4)
x T m and x Rn are respectively the transmitting and the receiving antenna positions, θ 0 refers to the signal Direction of Arrival (DoA) in azimutal plane considering the target at farfield distance from the radar system and k the wavenumber. The received signal for each frequency sample y in the setup depicted in Fig. 1 is thus given by the equation bellow:

y = H T s + ε (5)
where s and ε are respectively the transmitted signal vector and the Gaussian additive noise. An estimation of channel matrix H r can be retrieved computing the pseudo-inverse s + of the transmitted pseudo-orthogonal signals s:

H T r = ys + = H T ss + + ε × s + (6)
Considering the pseudo-orthogonality of the radiated waveforms, the following relation holds :

ss + ≈ [I] M×M (7)
leading to a perfect estimation of H in a noiseless scenario. MIMO imaging algorithms [START_REF] Bliss | Multiple-input multiple-output (MIMO) radar and imaging: degrees of freedom and resolution[END_REF]- [START_REF] Kpré | Experimental implementation of an ultra-wide band MIMO radar[END_REF] can then be applied to reconstruct the target image using the channel matrix H. Since the transmitted waveforms should be orthogonal for a simultaneous measurement of the channel matrix, each antenna should be addressed by an independent single waveform. As consequence, the number of transmitters increases with the number of transmit antennas needed.

In this context, the main contribution of this paper is introduction of a comparable technique achieved in the physical layer, allowing for a simplification of the MIMO RF architectures using a single transmitter to feed the M transmit antennas. The following section introduces this principle.

III. SINGLE INPUT MIMO RADAR PRINCIPLE

The technique presented here is based on a single input transmit array. The orthogonality of the transmitted waveforms is reached by the propagation in a 1 × M microwave cavity connected to the antennas. A generic setup is described in Fig. 2.

s M (f ) cM (f ) s m (f ) cm ( f ) s 1 (f ) c 1 ( f ) . . . . . .

x(f )

In p u t T r a n s m it te r 1 Fig. 2. Illustration of MIMO Radar transmitters feeding with 1 × M passive microwave cavity. Only one RF chain is needed in opposition with the conventional MIMO transmitter array.

The transmitted signals s(t) must be pseudo-orthogonal to estimate the channel matrix H. This requirement imposes a constraint on the component transfer functions which should be uncorrelated. This feature is the key factor of the single port MIMO radar. This is reachable with a metallic microwave oversized cavity which presents a high frequency diversity [START_REF] Carsenat | UWB Antennas Beamforming Using Passive Time-Reversal Device[END_REF]. In the depicted Fig. 2, the signal s m ( f ) radiated by the antenna m is :

s m ( f ) = x( f ) c m ( f ) (8)
where x( f ) is a signal with good auto-correlation properties generated at the input of the component, and c m ( f ) the component m th transfer function. For the sake of notation simplicity, Eq. ( 8) can be presented in a matrix form:

s = x c (9)
Therefore, the signals measured by the receiving antenna array considering a noise less scenario is:

y = H T s (10)
with H the previously defined channel matrix which accounts for useful informations required for the microwave imaging. An estimation of the channel matrix can be computed at each frequency by compensating the contribution of the generated waveform x and the component transfer function c in a single operation:

H T r = ys + (11)
where s + is the transmitted waveform pseudo-inverse. Numerous regularization techniques can be applied to solve this linear problem, such as Tikhonov regularization giving the following pseudo-inverse expression:

s + µ = (s † s + µI) -1 s † (12)
where µ denotes the regularization parameter which is used as a threshold to prevent any ill-conditioning problem [START_REF] Fromenteze | Waveform Coding for Passive Multiplexing: Application to Microwave Imaging[END_REF]. (.) † denotes the transposed-conjugate operator. The expression of the estimated channel matrix can be expressed by:

H T rµ = ys + µ = H T ss + µ ≈ H T R µ (13)
Considering an ideal equalization of the generated signal x, R µ can be considered as the correlation matrix of the transfer functions:

R µ =     c 1 × c + 1µ c 1 × c + 2µ . . . c 1 × c + mµ c 2 × c + 1µ c 2 × c + 2µ . . . c 2 × c + mµ . . . . . . . . . . . . c m × c + 1µ c m × c + 2µ . . . c m × c + mµ     (14)
Ideally, this matrix tends to an identity matrix, corresponding to a perfectly orthogonal transfer functions leading to a perfect estimation of the channel matrix H T rµ = H T . Such constraints are not practically achievable, limited by the frequency diversity of the chaotic component, but can be approached using a microwave cavity with regards to the operating wavelengths [START_REF] Kildal | Correlation and capacity of MIMO systems and mutual coupling, radiation efficiency, and diversity gain of their antennas: simulations and measurements in a reverberation chamber[END_REF].

Therefore, only a single wideband signal is generated at the input of a chaotic microwave cavity. The pseudo-orthogonality of radiated waveforms is ensured by the pre-characterized device transfer functions. The next section presents a numerical comparison between a conventional MIMO array and the proposed method.

IV. NUMERICAL RESULTS

A simulation of a radar experiment is presented to study the proposed method. The radar scenario consists of M = 4 transmit and N = 4 receiving antennas with inter-elements spacing respectively of dT = 1.4λ c and dR = 0.7λ c . Where λ c is the central wavelength in the 2.5Ghz to 3.5Ghz band. The inter element spacing are chosen in order to exploit MIMO radar virtual array [START_REF] Kpré | Experimental implementation of an ultra-wide band MIMO radar[END_REF]. A point-like target with an isotropic Radar Cross Section (RCS) is placed in the radar field of view (FOV) at 2m in range and 0.5m in cross-range which corresponds to an angle offset of 14 • in azimuth. The first simulation setup is a conventional MIMO radar in which the transmitted waveforms are modulated Q = 32 carriers and are computed via Eq.1.

In the second scenario, the chaotic cavity impulse responses are simulated using Rayleigh distributed magnitude modulated by an exponential decay with a time decay τ = 0.4µs which account for the quality factor of the cavity [START_REF] David | Hill Electromagnetic Theory of Reverberation Chambers[END_REF]. This simple mathematical model is used to reproduce the general behavior of an oversized cavity, leading to the following transfer function addressing the m th antenna :

c m ( f ) = F [r m (t) exp(-t/τ)] ( 15
)
where r m (t) denotes a random Gaussian noise sequence and F is the Fourier transform. Fig. 4 shows the comparison between the conventional MIMO radar and the proposed MIMO technique. Using the same bandwidth and the same post-processing technique. The results show that the proposed method can be an efficient alternative for the conventional MIMO architecture simplification. Despites a SNR reduction (4dB compared to the conventional MIMO radar), the compressed MIMO method can reach the same performances of the conventional with the main advantages of saving orthogonal waveforms generation and architecture simplification. The resolution of the beam at -3dB of the maximum is the same for both techniques (8.5 • ). In this simulation conditions the dynamic (main lobe to the higher sidelobe ration) is almost the same in azimuth and is about 16dB. This setting depends mainly of the correlation levels of the component transfer functions.

V. MEASUREMENT RESULTS AND DISCUSS

A measurement setup is set to match the simulation. An empty metallic microwave cavity with outer dimensions of 0.8×0.8×1m 3 is used to connect M = 4 transmitters. The N = 4 receivers are directly connected to an oscilloscope Agilent DSA90404A 20 GSa/s. The inter elements spacing is dT = 1.4×λ c for the transmitting and dR = 0.7×λ c for the receiving antennas. Fig. 5 shows the measurement setup. The metallic cavity impulse responses have been precharacterized as shown in Fig. 6. Good correlation properties are reached thanks to the frequency diversity features of the cavity.

An Arbitrary waveform generator Agilent M8190A 12 GSa/s is used to generate the single waveform (from Eq. 1) at the input of the cavity, and the oscilloscope is used to directly measure the echo signals. A watter bottle is set in front of the radar as shown in Fig. 7.

The backscattered signals are measured and post-processed to retrieve the bottle position. Fig. 6.b-c show the beamforming results which uphold the feasibility of the proposed method. With a single generated waveform, the MIMO channel matrix can be measured in one shot. The beam width at -3dB of the maximum is about 8 • which tends to the theoretical one. The dynamic is 12dB in the azimuth and the SNR of the image is about 18dB. This SNR is computed considering the real image and the background noise ratio.

VI. CONCLUSION

In this paper, a technique implemented for the simplification of MIMO radar architectures has been proposed. Using a passive microwave component whose transfer functions are quasi-orthogonal, MIMO radar transmit antennas can be independently addressed via a single generated waveform. This condition has been achieved with an oversized microwave cavity in regards to the operating wavelengths. The simulation results show that the proposed method can be applied to save orthogonal waveforms generation and to minimize the complexity of the conventional MIMO architecture while maintaining the same number of antennas. Future works will focus on the design of miniaturized microwave components that can provide good correlation properties. Herein the component is applied only at the transmit array, but it can also be used at the receiving array such that only a unique signal is measured at the reception leading to a 2 × MN ports passive component feeding M × N MIMO transmit/receiving antennas.

Fig. 1 .

 1 Fig. 1. Illustration of a conventional M × N MIMO radar. Each transmit antenna is fed by a transmitter. x T m and x Rn are respectively the transmit and receive antenna positions.

Fig. 3 .

 3 Fig. 3. MIMO Radar array. The virtual array is the convolution of the real transmit and receiving antennas aperture. The undermentioned numbers represent virtual elements redundancy.

Fig. 4 .

 4 Fig. 4. Comparison of the time domain beamforming in the conventional MIMO case and the proposed one.

Fig. 5 .

 5 Fig. 5. Single input MIMO radar setup. (a) UWB Tapered Vivaldi antennas array. The transmitters are connected to the metallic cavity behind and the receivers are connected to an oscilloscope. (b) Measuring and post processing bench.

Fig. 6 .

 6 Fig. 6. Metallic cavity general behavior. (a) 1 rst channel impulse response. (b) 1 rst channel autocorrelation r 11 and cross-correlation r 12 function between channel 1 and 2.

Fig. 7 .

 7 Fig. 7. Proposed MIMO radar measurement imaging results. (a) Radar scenario, a water bottle placed in front the Radar.(b) MIMO beamforming image. (c) Beam pattern at the target range (R = 1.8) plane.


Florent Jaillet
email: florent.jaillet@genesis.fr

Bruno Torrésani
email: bruno.torresani@cmi.univ-mrs.fr

Adaptive time-frequency representation for sound analysis and processing *

Introduction

The readability of a time-frequency representation generally depend crucially on a priori choices of some analysis parameters, which are often chosen heuristically. We describe here a set of methods for automatically selecting optimal parameters, exploiting sparsity requirements. Using Shannon or Renyi entropies for defining sparsity, we show that the "optimal" representation unfortunately often depends on the criterion. Therefore, we present several ways out for correcting such a shortcoming. In particular, we exploit the idea of "local optimization" in the time-frequency plane. Two different approaches are presented. In a supervised approach, the user can manually select the time-frequency domain in which window optimization is to be performed. In the unsupervised case, an iterative algorithm yields an automatic time-frequency segmentation, together with a corresponding inversion algorithm. As a by-product, it also provides multilayered representations for signals as in [START_REF] Daudet | Hybrid representations for audiophonic signal encoding[END_REF] [4], using reconstructions from a given window type only.

Sparsity measure

We introduce in this section the main tools we shall be using thereafter. For the sake of simplicity, we shall avoid here discretization issues, and describe our approach in terms of the short time Fourier transform (STFT), defined as follows [START_REF] Carmona | Practical Time-Frequency Analysis, volume 9 of Wavelet Analysis and its Applications[END_REF]. Given a finite energy window g (assuming that g = 1), associate with the finite energy signal x the function G g

x on the time-frequency plane

G g x (τ, ν) = x(t)g(t -τ ) e -2iπν(t-τ ) dt . (1)
As is well known, the STFT is invertible:

x(t) = G g x (τ, ν)g(t -τ )e 2iπν(t-τ ) dτ dν , (2)
and preserves energy

|G g x (τ, ν)| 2 dτ dν = |x(t)| 2 dt . (3)
Clearly, different windows g yield different time-frequency representations of the same signal. We are concerned here with the problem of selecting adaptively the "optimal" time-frequency representation for a given signal.

To do so, we have to introduce a criterion for optimality, and sparsity is a fairly natural choice: the idea is to find the representation which optimally concentrates the energy of the signal on a small number of time-frequency atoms. Unfortunately, there does not exist a unique way of measuring sparsity, and we shall see below that different sparsity measures yield different "optimal" representations. According to classical choices [START_REF] Wickerhauser | Adapted Wavelet Analysis from Theory to Software[END_REF], we shall limit ourselves to the family of Renyi entropies for defining sparsity. Consider a finite energy signal x, its short time Fourier transform G γ x with window γ, and introduce its

normalized spectrogram ρ γ (τ, ν) = |G γ x (τ, ν)| 2 / x 2 . For α ∈ (0, 1), the corresponding α-Renyi entropy reads R α (γ) = 1 1 -α log ρ γ (τ, ν) α dτ dν , (4)
and the Shannon entropy is obtained as the limiting case

S(γ) = lim α→1 R α (γ).
Given a family of window functions γ ∈ Γ, the "optimal" one will be defined as the one which minimizes the chosen entropy. We shall more specially focus on window families Γ consisting of dilates of a single window: γ s (t) = s -1/2 g(t/s), therefore looking for an optimal time-frequency representation from libraries of timefrequency atoms that have been considered by several authors in the literature [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF]. Although we have only presented here the discussion in the continuous case, the same argument goes through directly in the more practical situation where discrete Gabor expansions are considered instead of the STFT.

Supervised adaptation

Except for some specific cases, the size of the optimal window will depend highly on the entropy chosen to define the optimality. For example, using R α with a signal containing time-frequency atoms of different characteristics, smaller value of α turn out to yield smaller optimal windows. This shows that the notion of optimal STFT representation is not well defined for signals containing components with various time-frequency characteristics. But for a signal containing a well localized time-frequency component, explicit calculation on simple cases and numerical experiments show that the influence of the criterion on the optimal choice is higly reduced or even negligible. We thus introduce a first step in the adaptation to insure the localization of the components of the signal, by limiting the optimization to a user-defined region of the time-frequency plane. More precisely, for a given region Ω, using a given analysis window g, the corresponding signal x g,Ω is reconstructed with

x g,Ω (t) = Ω G g x (τ, ν) g(t -τ ) e i2πν(t-τ ) dτ dν . ( 5
)
The selection of the optimal window is then performed by minimizing the chosen entropy for this new signal x g,Ω .

Unsupervised adaptation

A natural extension of the latter approach consists in introducing a prior partitioning of the whole time-frequency plane and doing the adaptation in each subdomain. A simple example of such an approach can be built using only two windows, a "narrow" one g and a "wide" one h. This case is of special interest for the analysis and processing of audio signals. Indeed, these signals often contain two main class of components with different timefrequency characteristics. On one side, transient components are very localized in time and spread in frequency, and so best represented using a narrow window. On the other side, tonal components are slowly varying in time and well localized in frequency and so represented using a wide window. Therefore, as the automatic adaptation determines which components of the signal are best represented using a narrow or a wide window, it gives the opportunity to separate transient and tonal components and to represent each with suitable parameters. This algorithm could then for example be used as a precomputation to improve the quality of decompositions of the type described in [START_REF] Serra | A System for Sound Analysis/Transformation/Synthesis based on a Deterministic plus Stochastic Decomposition[END_REF].

To perform the unsupervised adaptation we use an iterative algorithm to control the reconstruction error. We first define a tiling of the time-frequency plane into rectangular "super-tiles" defined by

m,n = [m∆ τ , (m + 1)∆ τ ) × [n∆ ν , (n + 1)∆ ν ) , (6)
with ∆ τ > 0 and ∆ ν > 0 the time and frequency widths of the tiles. We introduce localized versions of the previously defined entropies. For a given α ∈ (0, 1), for a signal y and a window γ, we define

C m,n y (γ) = 1 1 -α log m,n |G γ y (τ, ν)| 2 E m,n γ (y) α dτ dν (7) with E m,n γ (y) = m,n |G γ y (τ, ν)| 2 dτ dν the energy in m,n .
The algorithm is initialized by defining r (0) = x. At each iteration, starting from k = 1, the set of tiles for which the window g is better than the window h

S k = {(m, n) ∈ Z 2 | C m,n r (k-1) (g) < C m,n r (k-1) (h)} (8)
is determined. From this, the corresponding partial reconstruction x

(k) g of the residual r (k-1) is obtained via

x (k) g = (m,n)∈S k m,n G g r (k-1) (τ, ν) g(t -τ ) e 2iπν(t-τ ) dτ dν, (9)
and from the complementary subset Z 2 \ S k the other partial reconstruction x

(k) h of r (k-1) reads x (k) h = (m,n)∈Z 2 \S k m,n G h r (k-1) (τ, ν) h(t -τ ) e 2iπν(t-τ ) dτ dν.
(10)

x (k) = x (k) g +x (k)
h is then an approximation of r (k-1) which is used to construct the new residual r (k) = r (k-1) -x (k) for the next iteration. Finally, we introduce the two layers

x g,K = K k=1 x (k) g , x h,K = K k=1 x (k) h .
(11)

x g,K contains components that are best represented with the narrow window, i.e. mostly the transient part of x, and x h,K contains components that are best represented with the wide window, i.e. mostly the tonal part of x.

The approximation of x at iteration K is x K = x g,K +

x h,K . Numerical tests show that the approximation error decreases (and become very small) as K grows.

An example of separation obtained with this algorithm is shown on figure 1.

Figure 1 :

 1 Figure 1: Example of separation on a glockenspiel sound : spectrogram of the original sound (top), spectrogram of the transient part computed with the narrow window (middle), spectrogram of the tonal part computed with the wide window (bottom).

* This work was supported in part by the European Union's Human Potential Programme, under contratct HPRN-CT-2002-00285 (HASSIP)