Vladislav A Yastrebov
email: vladislav.yastrebov@mines-paristech.fr

Wave-filtering properties of elastically-asymmetric architected materials

Keywords: numbers: 46.55.+d, 81.05.Zx, 46.40.Cd, 46.40.Ff

We propose a class of metamaterials possessing an arbitrary strong elastic asymmetry, which emerges from contact interaction introduced in their architecture. The material can be made stiff in tension and soft in compression or vice versa. This asymmetry results in different wave celerities of tensile and compressive components of elastic waves. The faster component overtakes the slower one and results in their dissipative annihilation through high-frequency energy cascades. We discuss the mechanism of this damping, efficient signal-absorbing assemblies, and relevant material architectures.

Mechanical behaviour of elastically asymmetric materials depends on the "directionality" of the strain tensor, i.e. on the signs of principal strains. For certain orientations σ(ε) = -σ(-ε), where σ is the Cauchy stress tensor and ε infinitesimal strain tensors. Therefore, the Young's, shear and bulk moduli of such materials can depend on the loading direction. In addition to this asymmetry, elastic anisotropy can come in play. Such materials are rather common in nature and include granular materials, soils, materials with internal flaws (e.g. cracked rocks, concrete) and others [1]. In soft matter, living cells, due to elastic asymmetry of fibrin and collagen, are able to self-adjust in response to external loads [START_REF] Wang | Proc. Natl. Acad. Sci. U.S.A[END_REF][START_REF] Janmey | [END_REF][4][5]. However, for most solid materials the elastic asymmetry is rather small and is hard to control. The theoretical studies of quasi-static mechanical behaviour of such materials can be found in [1,[6][7][8][9]. The elastic wave propagation in asymmetric media was studied in [10][11][12][13]. Vibrational analyses of materials with internal impacts, which render them elastically asymmetric, were conducted in [14][15][16].

First, we propose a class of architected materials with elastic asymmetry, which emerges from internal contacts established between parts of the elemental cell (see Fig. 1). This asymmetry is controllable and inversible, non-destructive and arbitrary-strong, contrary to a marginal asymmetry occuring in natural solid materials. The idea behind this novel architecture is that the non-adhesive contact can bear only compressive loads and opens in case of tensile ones. This non-smooth and amplitude-independent non-linearity renders the internal contact a good candidate for widening or enhancing novel and meta properties of architectured materials (see e.g. [17][18][19]). For example, if a zero-thickness cut is introduced in the material, as shown in Fig. 1(a), the resulting elastic modulus along OX direction depends on the sign of deformation: in case of tension, only thin ligaments bear the load and thus macroscopically the material behaves as a soft one; on the contrary, in compression, the contact in the cut is closed and can fully bear the load, thus resulting in a stiffer elastic behaviour. The elastic asymmetry can be controlled by the dimensions of the cut and by materials used in the central and peripheral zones. For a simple design presented here, a rough estimation of the ratio of Young's moduli in tension and compression can be given from simple geometrical consideration [20] by

E + /E -≈ (1 -l i /l)E 2 /E 1 ,
where l is the square-cell size, l i is the side-length of Π-shape cut, and E 1,2 are Young's moduli of the inner and outer materials: to amplify the asymmetry, materials can be chosen such that E 1 > E 2 . Thus, the resulting asymmetry can greatly overpass the asymmetry of existing materials, which are stiffer in compression than in tension, such as rocks and concrete. The opposite asymmetry also occurs in Nature in fibrous materials [21][22][23][24] and living cells [4]: this asymmetry is based on local buckling of fibers under compression. In contrast to this mechanism, our architecture uses contact interaction to achieve comparable property. A novel stiff-in-tension and soft-in-compression architec-ture is presented in Fig. 1(b): the contact is activated in tension and renders the material stiffer than in compression; in the latter case only thin ligaments bear the load as far as the gap remains open [see Fig. 1(b)]. The asymmetry of this material can be also enhanced by combination of stiff and soft materials in the architecture.

Such materials demonstrate unusual properties in dynamics both in vibration and wave propagation [13,15,25]. The latter presents the main topic of this paper. Namely, we study propagation of elastic waves through a one-dimensional assembly of symmetric and asymmetric materials, the focus is put on the damping properties. We show that the elastic asymmetry modifies the energy dissipation mechanisms, ensuring rapid damping for even low-frequency signals.

A one-dimensional wave equation for longitudinal waves propagating through a bi-modulus elastic material [10] with viscous dissipation of Kelvin-Voigt type (an elastic spring connected in parallel with a damper) [START_REF] Eringen | Mechanics of continua[END_REF] is of the form:

ρu ,tt = E(u ,x + α|u ,x |) ,x + µu ,xxt , -1 < α < 1, ( 1
)
where u is the horizontal displacement, lower indices after comma denote partial derivation • ,x and • ,t with respect to coordinate and time, respectively; ρ is the mass density (kg/m 3 ), E is the elastic modulus (Pa), and µ is the viscosity (Pa•s). The adimensional factor α determines the material asymmetry: the elastic modulus is equal to

E + = E(1 + α) and E -= E(1 -α) in tension (u ,x > 0)
and compression (u ,x < 0), respectively. Thus, the tensile and compressive components of elastic waves propagate at celerities c ± = E ± /ρ, respectively. Note that for a signal which has either purely compressive or tensile deformation, the wave propagation is governed by a linear Kelvin-Voigt model, however the superposition principle for deformation of different signs does not hold. The behaviour is independent of the signal amplitude since the non-linearity is localized in a single point on deformation curve (change of sign), which is centred at zero deformation. For one-dimensional systems, elastic asymmetry reduces to bi-modulus material model, thus hereinafter these two terms will be used interchangeably.

We consider propagation of elastic waves in structures shown in Fig 1(c,d) which include one or two sections with bi-modulus materials. Oscillations are induced at an excitation point in an elastically symmetric section of a bar by applying a harmonic force f = f 0 sin(ω 0 t). Forcing frequency ω 0 is chosen in the interval in which viscous effects are almost negligible, i.e. µω 0 /ρc ± 2 ≪ 1. The induced elastic weakly-dispersive waves propagate to the left and to the right: on the left they are absorbed by an absorbing layer. On the right they pass through a single section of length L of bi-modulus material [Configuration 1, see Fig. 1(c)], which without loss of generality can be considered to be stiff-in-tension and soft-in-compression, i.e. α 1 > 0 and

c - c -< c + c - c - c -< c + c - c - c -< c + c - c + < c -
E + 1 > E - 1 .
In addition, we consider Configuration 2 [see Fig. 1(d)], which has an extra section of a bi-modulus material with the opposite asymmetry, i.e. it is soft-in-tension and stiffin-compression (α 2 = -α 1 < 0, E + 2 < E - 2 ). The bimodulus section(s), are followed by a symmetric elastic bar with another absorbing layer on its extremity. The equations of motion are solved numerically using finite differences and Störmer-Verlet integrator, see Supplementary material [START_REF] Yastrebov | Supplementary material[END_REF]. The spatial discretization unit l reflects the size of the architectured elemental cell, thus the boundaries between symmetric and asymmetric sections pass between two springs separated by the concentrated mass. If λ 0 ≫ l (where λ 0 is a signal wavelength), such a homogenized model represents an effective medium for asymmetric architectured materials with internal contacts. The model is less reliable for short wavelengths λ/l ∼ 1, but can be viewed as a first order approximation.

Consider Configuration 1 [Fig. 2(a,d)] with a single wave u ,t (x, t) = -v 0 sin(ω 0 (x/c 0 -t))H(tx/c 0 )H(x/c 0 -t + 2π/ω 0 ) propagating from left to right if ω 0 > 0 at celerity c 0 , H(•) here denotes the Heaviside function, which is used to isolate a single wave [START_REF]We used a Gaussian smoothing of signal tails to avoid high frequency oscillations[END_REF]. In this wave the tensile component is followed by a compressive one. We assume, that the elastic modulus of the symmetric section is equal to the compressive modulus of the bi-modulus material E 0 = E -and c 0 = c -.

When this single wave enters the bi-modulus section, the tensile component accelerates abruptly and propagates at celerity c + > c -, whereas the compressive component continues to propagate at the same celerity c -and passes smoothly from symmetric to bi-modulus section, and further to symmetric section. In contrast, the tensile component is partly reflected back towards the emitter because of the elastic contrast: it occurs in entering and escaping the bi-modulus section. Within the bi-modulus section, because of the difference in wave celerity of the compressive and tensile components, at reentering into elastically symmetric zone, these two components are separated by ∆t = L/(c + -c -) in time and by ∆x = L/(1 -c -/c + ) in space. These properties can be used to construct a wave filter, which (1) due to reflection can partly attenuate the passage of either tensile of compressive wave components, and (2) due to contrast in celerities can separate the tensile and compressive components in space/time.

If the order of wave components is reversed [Fig. 2(b,e)], i.e. ω 0 < 0, then the compressive component precedes the tensile one, and the system dynamics is more complex. The slow leading component is overtaken by the faster tensile one and they start to interfere. Note that contrary to the purely symmetric case considered in [13], the ratio of amplitudes of the tensile and compressive components propagating in the bi-modulus section is given by A + /A -= 2/(1+c + /c -) (due to reflection), and their wavelengths relate as λ -/λ + = c -/c + . The overlap process leads to emergence of high frequency oscillations and accompanying viscous dissipation, which results in partial or almost complete annihilation of tensile and compressive wave components [see Fig. 2(e)]. However, these oscillations do not necessarily imply highfrequency alteration between tension and compression. In the first stage they simply superpose with the compressive wave [29]. Since in Kelvin-Voigt model the amplitude of a wave with a real wavenumber k decays in time t as exp(-µk 2 t/2ρ); in the limit of high wavenumber k > 2c ± /µ, the waves are overdamped and the harmonic part fully disappears. At later stages, the oscillations produced by the overlap of tensile and compressive components lead to frequent alteration of the deformation sign. As known [14,15], oscillators with piecewise smooth characteristics posses sub-harmonic resonances at higher frequencies, which ensure relatively high amplitude and enhanced dissipation.

It is straightforward to find a propagation distance L o needed for the tensile and compressive parts of initially harmonic signal to superpose completely. We introduce the following notations: in time of tensile and compressive components, which is equal to a half of the period ∆T = π/ω 0 in the symmetric section (where the oscillations are forced). Equating the length needed for faster and slower waves to travel to the same spatio-temporal point and requiring the full-period overlap, we obtain the overlap distance

E max = max{E + , E -}, E min = min{E + , E -
L o = 2πc max ω 0 ( √ γ -1) (2)
If the symmetric section is stiff c 0 = c max , then the overlap distance is given by

L o = λ 0 /( √ γ -1); if c 0 = c min then L o = λ 0 + λ 0 /( √ γ -1).
The high-frequency cascades emerging in annihilation of tensile and compressive waves present a powerful dissipative and filtering mechanism. To test its properties we analyse the energy passing through the bi-modulus section(s) as a function of its length and forcing frequency ω 0 . The injected energy is computed as the work of the imposed force E in = T 0 f (t)u ,t dt/2 [30], where T is the forcing time. The transmitted energy is computed right after the bi-modulus section(s) [see detector in Fig. 1(c,d)] as accumulated energy E out passing through a single elemental volume E out = T2 T1 [E k,t + E p,t ]dt, where E k,t , E p,t are the power of kinetic and potential energies, respectively. The time interval [T 1 , T 2 ] is chosen such to grasp all waves and most of reflections excited by the incident signal. Because of the absorbing layer on the far-right end of the structure, the energy integrated over time is equal to the energy passed from left to right within the given time interval.

The ratio of transmitted to incident signal energy log 10 (E out /E in ) is plotted in Fig. 3(a,b) as a function of the forcing frequency ω 0 and the length of the bi-modulus section L, for Configurations 1 and 2, respectively. The following parameters were used E/ρ = 1 (m 2 /s 2 ), elemental cell length is l = 1 (m) and area A = 1 (m 2 ), α = 0.3, µ/ρ = 0.01 (m 2 /s 2 ). For Configuration 1, the x/l = 0 x/l = 100

x/l = 350

x/l = 600

x/l = 800 ω r 2ω r 3ω r FIG. 4: Spectral density Φ(ω) at locations marked in Fig. 2(d,e). (a) Configuration 1, scenario of Fig. 2(b,e). Initial peak at the incident frequency ω0 is visible, as well as an emerging high frequency peak at sub-resonance frequency 2ωr. High-energy cascade with an exponent in the interval [0.9, 1] is highlighted. (b) Spectral density remains almost constant for Configuration 1, scenario of Fig. 2(a,d) for separating tensile and compressive wave components.

transmitted energy reduces significantly if the length of bi-modulus section is greater than the overlap length L o given in Eq. ( 2). For Configuration 2, the same qualitative effect is observed, however, it is greatly amplified [31] by the bi-modulus section of the same length L but with the opposite asymmetry, which is introduced right after the first bi-modulus section. After the first (partial) annihilation of tensile and compressive waves, the remnants of the tensile component precedes the compressive ones [see also Fig. 2(d)]. The second section with the inverse asymmetry serves to collide them again and dissipate their energy. Multiple combinations of antisymmetric sections can be used to ensure even more efficient damping as long as within every section tensile and compressive components overtake each other, i.e. the length of the section is greater than the overlap length L ≥ L o (see Eq. ( 2)). From Eq. ( 2) and the simulation data it follows that the absorption is very efficient for high and low frequency as long as ω 0 > 2πc max /L( √ γ -1). A similar damping shall occur for a random wave packet, which contains a roughly equal proportion of tensile and compressive components which "annihilate" via the same mechanism.

The presence of high-frequency energy cascades can be shown through spectral analysis. In Fig. 4 the evolution of the spectral density of deformation u ,x is shown for the cases depicted in Fig. 2(a,d andb,e). It is computed at several locations x along the bi-modulus section [marked in Fig. 2(d,e)] for the whole time history as Φ(ω, x) = û,x û * ,x , where û,x (ω, x) = ∞ 0 e -iωt u ,x (x, t)dt is the temporal Fourier transform, • * denotes the conjugate value, and ω is the angular frequency. We assume that at t < 0 the studied system is at rest. The system is forced mainly at frequency ω 0 = 0.03 Hz, resulting in a smooth peak in Φ for the signal entering the bi-modulus section at x = 0. Note that the resonance frequencies of the bi-modulus element is given by [16]), where the main resonance occurs for j = 1 and highfrequency sub-harmonic resonances occur at j > 1 and j ∈ N. In the signal spectrum [Fig. 4(a)] for the case of wave overlap [Fig. 2(b,e)], a second peak emerges at j = 2 sub-harmonic frequency of elemental cells. This peak, which grows with the propagation distance, presents a sink for the energy transmitted from low frequencies. The energy spectrum thus contains two peaks connected via a power law decay section with the exponent between -1.0 and -0.9 as shown in Fig. 4(c). Note that the energy decay of the main frequency ω 0 , as well as the rise and subsequent decay of sub-frequencies nω 0 , for n ∈ N, n > 1 is consistent with the findings presented in [13], where the main features of this dissipation were grasped. In contrast, when tensile and compressive components separate [Fig. 2(a,d)], the spectrum depicted in Fig. 4(b) does not present any particularities and the energy content remains rather stable.

ω r = 2j √ 1 -α 2 /( √ 1 -α + √ 1 + α) E 0 /ρl 2 (see,
In conclusion, we proposed new architectured materials in which the elastic asymmetry can be finely adjusted by combining internal contacts and components of different stiffness. Further, we studied the propagation of one-dimensional elastic waves in the resulting elastically asymmetric media. We found that the overlap of tensile and compressive wave components propagating at different celerities results in emergence of high-frequency energy cascades. It can be seen as tensile-compressive wave annihilation, which presents a novel and powerful signal damping mechanism. Moreover, it is enhanced drastically if a section with one asymmetry is followed by another one with the opposite asymmetry. The efficient wave damping is ensured if the bi-modulus section is chosen longer than the wave-overlap length [START_REF] Wang | Proc. Natl. Acad. Sci. U.S.A[END_REF]. The ratio of the overlap length to the incident wavelength scales as L o /λ 0 ∼ 1/( √ γ -1). The key advantange of the proposed architected materials consists in fact that the elastic asymmetry γ can be adjusted to be high, which would enable to keep the damping device relatively small compared to the characteristic incident wavelength.

The demonstrated efficient damping mechanism can be used in shock absorbing systems, and, potentially, in seismic protection from surface waves. The one-dimensional model presented here can be extended to two and three dimensional cases, where the compressive/tensile elastic asymmetry should be enhanced with shear asymmetry and complemented by elastic anisotropy.

FIG. 1 :

 1 FIG. 1: (a,b) Examples of design of elemental architected cells with asymmetric elastic properties: the effective stiffness depends on the sign of horizontal strain component (compression/tension): architecture (a) is soft-in-tension and stiff-incompression, architecture (b) is stiff-in-compression and softin-tension. The associated deformation curves are also shown. Configurations 1 and 2 depicted in (c) and (d), respectively are used to study 1D wave propagation in asymmetric media.

FIG. 2 :

 2 FIG. 2: Spatio-temporal wave tracing is presented in (a) for a tensile component (dashed line) followed by a compressive one (solid line), (b) inverse incident signal, (c) inverse incident signal passing through a configuration with sections of opposite asymmetries; (a,b) and (c) correspond to Configuration 1 and 2, respectively [see Fig. 1(c,d)]. Simulation results for α = 0.3, µ/rho = 0.01 m 2 /s 2 , E/ρ = 1 m/s: (d) and (e) represent spatio-temporal deformation map (u,x) and correspond to the diagram (a) and (b), respectively; x/l = 0 separates symmetric and asymmetric sections.

FIG. 3 :

 3 FIG. 3: The decimal logarithm of the ratio of transmitted to incident energy is plotted for different forcing frequencies ω0 and different lengths of bi-modulus section L: (a) Configuration 1 [see Figs. 1(c),2(d)], (b) Configuration 2 [see Figs. 1(d),2(e)]. Dashed curves mark transmission iso-levels, dash-dotted curves represent Eq. (2).

(a) (b) (c) coordinate, x coordinate, x coordinate, x time, t time, t time, t(d)

Acknowledgement. The author is grateful to Lev Truskinovsky for valuable discussions, to Samuel Forest for his encouragement, and to Arsen Subashiev for helpful remarks to the text.


Dario Magliacano

Massimo Viscardi

Ignazio Dimino
email: i.dimino@cira.it

Antonio Concilio

ACTIVE VIBRATION CONTROL BY PIEZOCERAMIC ACTUA-TORS OF A CAR FLOOR PANEL

Active vibration control was successfully tested in many automotive and aeronautical applications. The main purpose of this work is the implementation and experimental investigation of the applicability and efficiency of an active vibration control (AVC) concept for vibration reduction in an automobile passenger compartment. In the contest of the presented activity, the floor panel of a medium-class test car has been studied as test case for an AVC implementation based on the use of piezoelectric actuators. A preliminary numerical activity has been performed for the correct evaluation of the dynamic behaviour of the system under the specific external disturbance spectra; based upon this information, a procedure for the correct positioning of sensors and actuators devices has been implemented. The following experimental activity has then been dedicated to the implementation of both SISO (Single Input Single Output) and SIMO (Single Input Multiple Output) feedforward control applications, under simulated conditions of external disturbance field. The principal results that will be discussed within the paper show a good ability of the actuation system to manage the control signals in terms of vibrational energy; this circumstance, candidate piezoelectric patches as adequate actuation devices. Under the specific control point of view, it appears evident that SISO results are appreciable only at error sensor, while uncontrollable conditions emerge at other locations. This limitation may be partially overcome with the SIMO architecture, but a review of the main limitations concerning this architecture are also presented.

Introduction 1.1 Active control applications

There is always a stronger and more insistent research, by producers, of new technical solutions and control methods for noise and vibration level reductions, with better performances and more flexible operational capabilities. Areas of application that now show an irreversible trend towards the adoption of active sound and/or vibration control systems are:

• Transport vehicles (aircrafts, cars, trains, ships), concerning the improvement of interior comfort and outdoor pollution; • Civil engineering structures (skyscrapers, towers, workshops, bridges, industrial complexes) regarding earthquake protection, the reduction of the effects due to noise and vibrations induced by road and rail traffic, by the wind, and by the presence inside and/or outside of large rotary machines, such as turbines, pumps or compressors; • Household appliances (refrigerators, washing machines and dishwashers).

Sources of interior noise and vibrations in automotive

With reference to the first class of systems and specifically, to automotive vehicle, it can be stated that a variety of independent noise sources impact on passengers' perceived noise in vehicles. These sources may be transmitted via a variety of structural paths, and then radiated acoustically into the cabin (structural-borne noise) or acoustically generated and propagated by airborne paths (air-borne noise). The major contributors may be thus distinguished in:

• Tire rolling noise: the interaction between tires and ground generates broadband vibrations that are transmitted to the structure. This phenomenon is more remarkable at low temperatures, due to the stiffening of the tire. Working on suspensions could severely affect performances, while the use of a structural shock absorber could heavily increase the weight; so, here is the need to develop an active control system. • Engine noise: mechanic and chemical processes inside the engine generate vibrations and therefore noise, which is transmitted to the cockpit. In addition, new low-power engines need a design that causes greatest noise at equal thrust. • Aerodynamic noise: the "noise vortices" are caused by the interaction between a fluid (air) and a non-aerodynamic body (car). The frequency of detachment of these Von Karman vortices generates unstable forces that produce both broadband and tonal noise. This noise is transmitted to the structure, that vibrates and transfers it into the interior.

Very often, the attenuation of these noise sources require a substantial increment of the mass and stiffness of structural elements present along the noise transmission path, with a consequent deterioration of other vehicle parameters (consumption, performances, drivability, ...); also because, these technologies are not able to adapt their performance to the specific condition and are so designed for the worst condition to be faced with. The engineering challenge of the lightweight structure need, by this point of view, new solution that may be partially represented by those recognized as "active technologies".

Active control with feedforward strategy 2.1 Introduction

The reduction of interior noise through active devices has been a research topic for more than 50 years. Rubber and hydro mounts are the standard tool to isolate the engine and the transmission from the chassis. Rubber isolators work well (in terms of isolation) when the rubber exhibits low stiffness and little internal damping. Little damping, however, leads to a large resonance peak which can manifest itself in excessive engine movements when this resonance is excited (front end shake). These movements must be avoided in the tight engine compartments of today's cars. A low stiffness, while also giving good isolation, leads to a large static engine displacement and to a low resonance frequency (which would adversely affect the vehicle comfort and might coincide with resonance frequencies of the suspension system). Classical mount (or suspension) design therefore tries to achieve a compromise between the conflicting requirements of acceptable damping and good isolation. It is clear that this, as well as other passive vibration control measures, is a trade-off design method in which the properties of the structure must be weighted between performance and comfort. An attractive alternative that overcomes the limitations of the purely passive approach is the use of active noise and vibration control techniques (ANC/AVC). The basic idea of ANC and AVC is to superimpose the unwanted noise or vibration signals with a cancelling signal of exactly the same magnitude and a phase difference of 180° (i.e. the ''anti-noise" principle of Lueg (Lueg, 1933)) [START_REF] Svaricek | Automotive Applications of Active Vibration Control[END_REF]. A feedforward control strategy is possible only if there is a preliminary information about the excitation of the primary source to be controlled; if that is not possible, feedback control strategies are preferred instead. Typically, the primary excitation cannot be directly measured and a reference sensor is used to detect the primary wave. The signal, from the reference sensor, is in turn transmitted to an electronic controller, which then controls the secondary excitation (control actuator); the main purpose is then to design a controller capable to minimize the output of one or many error sensors. More recently, the active controllers are implemented in digital signal processor (DSP) boards. This gives rise to the need for accounting inevitable delays associated with the controller, due to both the DSP processing time and the presence of amplifiers in the feedforward control scheme. The feedforward control of random noise, therefore, is effective only if the delay in the path from the reference sensor to the error (via the controller) is less than the time that the elastic wave takes to travel the same distance through the structure; this is generally referred to as "causality constraint". Conversely, this effect is largely manageable when the incident wave is of deterministic nature (e.g., harmonics). The current paper considers the reduction of car rear under-seat panel-induced vibrations. A pure feedforward strategy (with no adaptation), using piezoelectric patches as control actuators is implemented and experimented on the test article. The criterion of least squares (optimal control) is applied to estimate the performance of feedforward control systems using different local control schemes with one or more secondary sources and one or more error sensors.

Optimal control

Assuming a linear behaviour of all the parts of the system (structural, electrical and piezoelectric) and all the signals are in the respective stationary states, the frequency-domain complex output of the error sensor (y) is the complex output overlapping due to only the primary excitation (yd) and to exclusively the secondary one (yc). The error signal can then be expressed in terms of complex concentrated force amplitude (vd) and complex amplitude of the voltage applied to the piezoelectric (vc): 𝑦 = 𝑦 𝑑 + 𝑦 𝑐 = 𝐻 𝑑 𝑣 𝑑 + 𝐻 𝑐 𝑣 𝑐 In the previous equation, Hd represents the so-called primary path, i.e. the FRF from the primary source to the error sensor, while Hc is the secondary path, i.e. the FRF between the secondary source and the error sensor. It therefore assumes that the signals relating to primary and secondary sources are perfectly consistent, that is there are no other unrelated signals (such as noise); therefore, the acceleration signal can be cancelled at any frequency by introducing a secondary inputs equal to:

𝑣 𝑐 𝑚𝑖𝑛 = - 𝐻 𝑑 𝐻 𝑐 𝑣 𝑑
Obviously, the FRF of the secondary path should not be null to prevent an optimal control with divergent secondary input; this singular condition is avoided by introducing damping into the system. Generally, the control system consists of multiple secondary sources and multiple error sensors (MIMO system); in this case, the equation can be rewritten as: {𝑦} = [𝐻 𝑑 ]{𝑣 𝑑 } + [𝐻 𝑐 ]{𝑣 𝑐 } Note that the (i, j)-th element of the [Hc] matrix is the FRF between the i-th error sensor and the j-th secondary source. If the number of error sensors and the number of secondary sources are equal, the optimal vector of the secondary inputs (i.e., the vector that eliminates the output of all error sensors) is given by:

{𝑣 𝑐 𝑚𝑖𝑛 } = -[𝐻 𝑐 -1 ][𝐻 𝑑 ]{𝑣 𝑑 }
In many cases, however, the number of error sensors and that of secondary sources are not equal; in such cases, the optimal secondary input can be calculated by minimizing the error through the application of the least squares criterion, in the form: 𝐽 = {𝑦 𝐻 }[𝑊 𝑦 ]{𝑦} + {𝑣 𝑐 𝐻 }[𝑊 𝑐 ]{𝑣 𝑐 } Where [Wy] is the weight matrix of the errors and [Wc] is the gain weight matrix; through the first, some of the error signals may have priority over others. The error criterion can be written in the socalled "Hermitian quadratic form": 𝐽 = {𝑣 𝑐 𝐻 }[𝐴]{𝑣 𝑐 } + {𝑣 𝑐 𝐻 }{𝑏} + {𝑏 𝐻 }{𝑣 𝑐 } + 𝑐 Where single terms are given by:

[𝐴] = [𝐻 𝑐 𝐻 ][𝑊 𝑦 ][𝐻 𝑐 ] + [𝑊 𝑐 ] {𝑏} = [𝐻 𝑐 𝐻 ][𝑊 𝑦 ][𝐻 𝑑 ]{𝑣 𝑑 } 𝑐 = {𝑣 𝑑 𝐻 }[𝐻 𝑑 𝐻 ][𝑊 𝑦 ][𝐻 𝑑 ]{𝑣 𝑑 }
The cost function J has a unique minimum if [A] is not singular and this is guaranteed if [A] is positive defined; to verify this condition, the weight matrices [W y ] and [W c ] must be positive defined. Downstream of this, the control inputs that minimize the error criterion are: {𝑣 𝑐 𝑚𝑖𝑛 } = -[𝐴 -1 ]{𝑏} While the minimum of the error criterion is:

𝐽 𝑚𝑖𝑛 = 𝑐 -{𝑏 𝐻 }[𝐴 -1 ]{𝑏}

System dynamic characterization and actuator's location

The test article of the present activity, as previously stated, is represented by the rear floor panel of a car; the squared box in figure 2 give evidence of the specific area. As a first activity, a numerical modal analysis of this component has been carried out from 0 Hz to 300 Hz using the software MSC Nastran; it revealed the presence of six normal modes of vibration within that range, whose relative frequencies are reported in next Table 1. It is important to underline that the numerical model include the mass and stiffness effect of the rest of the car body floor that has been condensed during a preliminary model refinement. Relative strain distribution contour plots are presented in the next pictures (figures 3-8).

Experimental setup

The experimental activities have been carried out on the rear car floor panel of the real car at the C.I.R.A. facilities. The primary excitation device is represented by a shaker (Figure 10), suspended with a hook and interfaced to the trunk of the car. Realistically, the primary excitation should originate from an area close to the suspensions; however, because of the difficulty of access to zone, downstream of the detection of a good transmissibility between the area of the trunk and that of the rear floor panel, it was opted for this solution.

Pure sine signals where used as excitation spectra during the first part of the activity; the main results herein presented refers to the 300 Hz pure tone (this frequency does not coincide with any normal mode and has been focused as significant frequency related to real operative conditions).

The error and monitoring sensor set-up consisted of 44 accelerometers installed on the rear floor panel (figure 9), even if has been noted that this acquisition and control grid has not been used, in its entirety, for the purpose of the tests discussed here, but it has been prepared as necessary for additional tests not covered in this paper. 5 locations were then selected for the piezoelectric patches (control sources), on the basis of specific strain contours at the selected control frequencies.

Control system consisted of a dSpace multichannel control board where two feedforward control algorithms have been implemented.

Experimental results

As already described in the introduction chapter, the present study was mainly targeted at the verification of the piezoelectric patches ability to introduce inside the system the adequate amount of energy to control the primary disturbance field preliminary evaluated on the basis of road test. To stress this verification a pair of locations were very far each other; in the specific, it was decided to use the piezoelectric "3" (circled in blue in Figure 9) as sensor and the piezoelectric "5" (circled in red in Figure 9) as control device.

To verify the ability of control signal to counteract the primary one the phase of the corresponding actuating signal has been calculated with the formula: control phase (in degrees) = 180 -(B2 -B1) where:

 B1 is the phase of the excitation devicesensor FRF;  B2 is the phase of the control devicesensor FRF.

Under this condition, it has been verified that a control voltage as small as 35 V was able to generate a secondary control field causing a strong reduction of overall vibration at the error sensor location. This circumstance represented an important preliminary results, because such a voltage level does not imply the use of big and expensive power supply but can be easily managed with light and cheap logic technology.

SISO application

For the first experimental application discussed further below, after the preliminary tests, a SISO architecture was implemented. During these test, different pair of sensors and actuators were tested. All the tests were performed under the application of a sinusoidal excitation signal, and as a general consideration it can be stated that, in the best condition, the controlled response was detected equal to 21 % of the not controlled one (Figure 12); in other words, the obtained reduction of the response equal has been always greater that 13.6 dB. This result appears satisfactory, but also expected, because under a theoretical point of view, the control system is always able to manage the right phase of the control signal to perfectly be in opposition of phase at the error sensor. Unfortunately, this circumstance may generate a possible increment of the vibration level at other location not managed within the control algorithm; see, in example what has been measured under this configuration where piezoelectric "3" (circled in blue in Figure 9) as error sensor, the piezoelectric "4" (circled in green in Figure 9) as monitoring sensor and the piezoelectric "2" (circled in yellow in Figure 9) as control device. Under the above conditions, it has been checked an amplification of +22.8 dB of the vibrational level in correspondence of the monitoring sensor; but also this circumstance was expected because both the optimum phase or amplitude parameters referred to the error sensor may have a complete opposite effect at other locations (in this case, in example, the increment is more significant at monitoring station than the reduction at the error location, due to the proximity to the actuation source).

Figure 13: Effect of the SISO feedforward control at error and monitoring sensors.

SIMO application

The second application, in addition to the above, provides the presence of two error sensors (both sensor 3 and 4 have been introduced in the control low); in this circumstance, a reduced level of the controlled response may be measured at both locations, but relative value (approx. -2.5 dB in both sensors) is much reduced if compared at SISO configuration. Also this behaviour was expected because the searched minimum vibration level at both sensors requires a compromise of the optimal amplitude and phase of the control signal.

Conclusions and future planned activities

The present study gathered essential data and parameters for the design of noise and vibration active control architectures in the automotive sector. Preliminarily numerical and experimental activities have shown the opportunity to use piezoelectric devices both as sensor and actuators because they have shown a good sensibility when used as sensors and an adequate transmitted energy when used as actuators. SISO and SIMO test herein presented have to be considered as preliminary approach to validate the physical architecture; it is evident that a MIMO feedforward approach or a feedback strategy will be realistic implemented on the basis of availability of a correlated reference signal or not.

Figure 1 :

 1 Figure1: Feedforward control system to control the propagation of a wave within an elastic structure.[START_REF] Nijhuis | Analysis tools for the design of active structural acoustic control systems[END_REF]

Figure 2 :

 2 Figure 2: Car floor panel.

Figure 3 :

 3 Figure 3: Strainsmode 1.Figure4: Strainsmode 2.

Figure 4 :

 4 Figure 3: Strainsmode 1.Figure4: Strainsmode 2.

Figure 5 :

 5 Figure 5: Strainsmode 3.Figure6: Strainsmode 4.

Figure 6 :

 6 Figure 5: Strainsmode 3.Figure6: Strainsmode 4.

Figure 7 :

 7 Figure 7: Strainsmode 5.Figure8: Strainsmode 6.

Figure 8 :

 8 Figure 7: Strainsmode 5.Figure8: Strainsmode 6.

Figure 9 :

 9 Figure 9: Experimental setup.Figure10: Excitation device.

Figure 11 :

 11 Figure 11: Signal Fourier Transform.

Figure 12 :

 12 Figure 12: Effect of the feedforward control (left) and zoom (right).

Figure 14 :

 14 Figure 14: Effect of the feedforward control at error sensors 3 and 4.

Table 1 :

 1 Numerical proper modes of vibration.

	Mode	1	2	3	4	5	6
	Frequency [Hz]	152.79	204.56	236.21	259.39	259.49	283.11


Matteo Conforti
email: matteo.conforti@ing.unibs.it

Fabio Baronio

Extreme high-intensity and ultra-broadband interactions in anisotropic β-BaB 2 O 4 (BBO) crystals

Keywords:

We derive unidirectional pulse propagation equations to describe extreme high-intensity and ultra-broadband optical interactions in uniaxial crystals, showing both second-and third-order nonlinear optical susceptivities. We focus our attention on the anisotropic nature of the quadratic and cubic nonlinear response of β-BaB 2 O 4 (BBO) crystals. Two nonlinearly coupled first order (in the propagation coordinate) equations describe the dynamics and interactions of the ordinary and extraordinary field polarizations, and are valid for arbitrarily wide pulse bandwidth. We exploit this model to predict harmonic and supercontinuum generation in BBO crystals under strong and competing influence of quadratic and cubic susceptivities.

I. INTRODUCTION

In recent years there has been a great deal of interest in research on second-harmonic (SHG) [START_REF] Cittadini | I like fuck Fabio[END_REF], high-order harmonic (HHG) [START_REF] Krausz | Attosecond physics[END_REF], and supercontinuum (SC) generation [START_REF] Dudley | Supercontinuum generation in photonic crystal fiber[END_REF] in nonlinear optical media for such diverse applications as frequency metrology, few-cycle pulse generation, spectroscopy, biological and medical analyses.

The SHG of super-strong ultrashort (tens of femtoseconds) laser pulses, using the χ (2) nonlinearities in optical crystals, is a very important task, because the process can be used not only for wavelength conversion, but for significant improvement of temporal intensity contrast ratio and pulse shortening. SHG is especially important for Ti:sapphire laser facilities operating at 800 nm [START_REF] Aoyama | 0.85-PW, 33-fs Ti:sapphire laser[END_REF] and optical parametric amplifiers at 910 nm [START_REF] Lozhkarev | Compact 0.56 petawatt laser system based on optical parametric chirped pulse amplification in KDP crystals[END_REF].

SC generation has been performed conventionally using the χ (3) nonlinearities in optical fibers. Due to the high nonlinearity and engineerable dispersion available in fibers, spectra spanning multiple octaves can be achieved [START_REF] Farrell | Octave-spanning super-continuum from a silica photonic crystal fiber pumped by a 386 MHz Yb:fiber laser[END_REF][START_REF] Fang | Multiwatt octave-spanning supercontinuum generation in multicore photonic-crystal fiber[END_REF]. However, reaching the mid-infrared spectral region with χ (3) -based SC sources is challenging [START_REF] Price | Mid-IR supercontin-uum generation from nonsilica microstructured optical fibers[END_REF]. A promising alternative approach consists on the exploitation of the χ (2) nonlinearities of optical crystals for SC generation [START_REF] Conforti | Nonlinear envelope equation for broadband optical pulses in quadratic media[END_REF][START_REF] Phillips | Supercontinuum generation in quasi-phase-matched LiNbO3 waveguide pumped by a Tmdoped fiber laser system[END_REF]. SC interactions can readily be achieved in birefringent or quasi-phase matched (QPM) crystals [START_REF] Conforti | Ultra-broadband optical phenomena in quadratic nonlinear media[END_REF][START_REF] Phillips | Supercontinuum generation in quasi-phasematched waveguides[END_REF], with high-intensity light pulse excitation. Quadratic SC generation, difference frequency generation and optical parametric generation are currently active areas of research [START_REF] Zhou | Ultrafast and Octave-Spanning Optical Nonlinearities from Strongly Phase-Mismatched Quadratic Interactions[END_REF][START_REF] Levenius | Multistep quadratic cascading in broadband optical parametric generation[END_REF].

Nowadays, technological advances in ultrafast optics have permitted the generation of ultraintense light pulses comprising merely a few field oscillation cycles. Peak intensities approaches 10 15 W/cm 2 [START_REF] Sung | 0.1 Hz 1.0 PW Ti:sapphire laser[END_REF], opening the study of an entirely new realm of nonlinear interactions in solid materials.

Beta-Barium-Borate (β-BaB 2 O 4 , BBO) is a very popular crystal, among all solid-state optical materials: BBO has a high damage threshold, low dispersion and χ (2) nonlinearities of few pm/V allowing for efficient quadratic frequency conversion interactions [START_REF] Nikogosyan | Nonlinear Optical Crystals: A Complete Survey[END_REF].

In this work, we explore the use of BBO crystals in extreme optical regimes, where dispersion effects and cubic nonlinearities play an essential role. In particular, we derive a comprehensive model to describe the propagation of extreme high-intensity and ultrabroadband optical pulses in BBO crystals. This model provides a powerful tool due to its generality and simplicity, and can be easily solved with a modest computational effort.

The paper is organized as follows. In Section 2, we recall the derivation of the master equations in uniaxial media, discussing the validity of the model. We consider both the second-and third-order nonlinear contributions, and their angular dependences. We take into account all possible second-and third-order interactions, including ones typically nonphase-matchable. In Section 3, we present some numerical examples of second harmonic generation and supercontinuum generation in BBO crystals, showing the key role of cubic susceptivity. Eventually we draw our conclusions in Section 4.

II. DERIVATION OF THE PROPAGATION EQUATIONS

In this section we review and extend the derivation of the unidirectional nonlinear vector field equations reported in [START_REF] Conforti | Modeling of ultrabroadband and single-cycle phenomena in anisotropic quadratic crystals[END_REF] (also called Forward Maxwell Equations, FME [START_REF] Housakou | Supercontinuum generation of higher-order solitons by fission in photonic crystal fibers[END_REF], or Unidirectional Pulse Propagation Equation, UPPE [START_REF] Kolesik | Unidirectional optical pulse propagation equation[END_REF]), describing the propagation of the ordinary and extraordinary polarizations of the electric field in uniaxial crystals with both χ (2) and χ (3) nonlinearities.

We start from Maxwell equations written in MKS units, in the reference frame x y z

∇ × E = - ∂B ∂t (1)
∇ × H = ∂D ∂t (2)
B = µ 0 H (3)
D = D L + P N L (4)
where D L and P N L account for the linear and nonlinear response of the medium, respectively. The components of the linear displacement vector for a dispersive anisotropic medium reads (assuming summation over repeated indxes)

D L,j = ε 0 ∞ -∞ ε jk (t -t )E k (t )dt . (5)
In the reference frame of the principal axes of a uniaxial crystal, the dielectric permittivity tensor is the diagonal matrix ε = diag(ε o , ε o , ε e ), where ε o , ε e are the ordinary and extraordinary relative dielectric permittivity, respectively. The reference frame of the principal axes of the crystal (x y z ) is not convenient for the derivation of the propagation equations.

We introduce a reference frame xyz that is rotated by (θ, φ) with respect to crystal axes.

Namely, θ is the angle between the propagation vector (parallel to z) and the crystalline z axis (the crystal optical axis), and φ is the azimuthal angle between the propagation vector and the x z crystalline plane. The two reference frame are linked by the orthogonal rotation matrix A:

A =       cos φ cos θ sin φ cos θ -sin θ -sin φ cos φ 0 sin θ cos φ sin φ sin θ cos θ       . ( 6
)
The dielectric permittivity tensor in the xyz frame is no longer diagonal, and it can be written as

ε = Aε A T =       ε o cos 2 θ + ε e sin 2 θ 0 (ε o -ε e ) cos θ sin θ 0 ε o 0 (ε o -ε e ) cos θ sin θ 0 ε o sin 2 θ + ε e cos 2 θ       . (7)
In the reference frame xyz, it is possible to decompose the electromagnetic field into two linear and orthogonal polarizations of D, both transverse to the propagation direction

z [20]: D = (0, D y , 0) T + (D x , 0, 0) T .
We assume the propagation of plane waves, so the electric field and displacement vectors depend upon the z coordinate (and time) only. It is worth noting that this decomposition is rigorous for linear propagation only, since the nonlinearity can rotate locally the polarization. However it is reasonable to consider the nonlinearity as a perturbative term whose effect is to couple the orthogonal polarized field vector components during propagation. If we neglect dispersion and nonlinearity, just for the moment, the electric field vector can be straightforwardly computed as:

E = ε -1 0 ε -1 D = ε -1 0       cos 2 θ εo + sin 2 θ εe D x ε -1 o D y εe-εo εeεo cos θ sin θD x       (8)
By eliminating the magnetic field from Maxwell equations we obtain the vector wave equation:

∇ × ∇ × E - 1 ε 0 c 2 ∂ 2 D L ∂t 2 = 1 ε 0 c 2 ∂ 2 P N L ∂t 2 (9)
Note that obviously ∇ • D = 0, but ∇ • E = 0. By writing (9) in components we obtain

∂ 2 E x ∂z 2 - 1 ε 0 c 2 ∂ 2 D L,x ∂t 2 = 1 ε 0 c 2 ∂ 2 P N L,x ∂t 2 (10) ∂ 2 E y ∂z 2 - 1 ε 0 c 2 ∂ 2 D L,y ∂t 2 = 1 ε 0 c 2 ∂ 2 P N L,y ∂t 2 (11) 0 = 1 ε 0 c 2 ∂ 2 P N L,z ∂t 2 (12)
The last equation witnesses the fact that the decomposition into two independent orthogonal polarizations is rigorous only in the linear case. We neglect P N L,z , in the reasonable hypothesis of small nonlinearity.

Exploiting the relation ( 5) we obtain:

∂ 2 E m (z, t) ∂z 2 - 1 c 2 ∂ 2 ∂t 2 +∞ -∞ E m (z, t )ε m (t -t )dt = 1 ε 0 c 2 ∂ 2 ∂t 2 P N L,m (z, t) , m = x, y (13)
where we have defined

ε x = cos 2 θ ε o + sin 2 θ ε e -1 (14)
ε y = ε o (15)
We thus have obtained the propagation equations for an ordinary polarized wave E y and an extraordinary polarized wave E x .

By defining the Fourier transform F[E](ω) = Ê(ω) = +∞ -∞ E(t)e -iωt dt, we can write (13) in the frequency domain:

∂ 2 Êm (z, ω) ∂z 2 + ω 2 c 2 εm (ω) Êm (z, ω) = - ω 2 ε 0 c 2 PNL,m (z, ω), ( 16
)
where c is the velocity of light in vacuum, ε 0 is the vacuum dielectric permittivity, εm (ω) = 1 + χm (ω), χm (ω) is the linear electric susceptibility and k m (ω) = (ω/c) εm (ω) is the propagation wavenumber.

We now proceed to obtain, from the second order vector wave equation ( 16), an equation, first order in the propagation coordinate z, describing electromagnetic fields propagating in the forward direction only. Several techniques have been proposed in literature in order to achieve a pulse propagation equation with minimal assumptions [START_REF] Housakou | Supercontinuum generation of higher-order solitons by fission in photonic crystal fibers[END_REF][START_REF] Kolesik | Unidirectional optical pulse propagation equation[END_REF][START_REF] Brabec | Intense few-cycle laser fields: Frontiers of nonlinear optics[END_REF][START_REF] Kolesik | Nonlinear optical pulse propagation simulation: From Maxwell's to unidirectional equations[END_REF][START_REF] Genty | Nonlinear envelope equation modeling of sub-cycle dynamics and harmonic generation in nonlinear waveguides[END_REF][START_REF] Kinsler | Theory of directional pulse propagation[END_REF][START_REF] Kinsler | Optical pulse propagation with minimal approximations[END_REF][START_REF] Kumar | Ultrashort pulse propagation in a cubic medium including the Raman effect[END_REF]. The interested reader can find in [START_REF] Kinsler | Optical pulse propagation with minimal approximations[END_REF][START_REF] Kolesik | Theory and simulation of ultrafast intense pulse propagation in extended media[END_REF] an exhaustive discussion on the different derivation styles. Here we decided to follow the approach outlined in the review paper [START_REF] Kolesik | Theory and simulation of ultrafast intense pulse propagation in extended media[END_REF], that combines minimal assumptions and straightforward derivation.

We write the electric field components in spectral domain as the sum of a forward (F ) and a backward (B) propagating part, that with our definition of the Fourier transform reads:

Êm (z, ω) = Fm (z, ω)e -ikm(ω)z + Bm (z, ω)e ikm(ω)z . ( 17
)
By plugging Ansatz ( 17) into ( 16), we get:

∂ 2 Fm ∂z 2 -2ik m (ω) ∂ Fm ∂z e -ikm(ω)z + ∂ 2 Bm ∂z 2 + 2ik m (ω) ∂ Bm ∂z e ikm(ω)z = - ω 2 ε 0 c 2 PNL,m ,
that can be rewritten as:

∂ ∂z ∂ Fm ∂z e -ikm(ω)z + ∂ Bm ∂z e ikm(ω)z - -ik m (ω) ∂ Fm ∂z e -ikm(ω)z - ∂ Bm ∂z e ikm(ω)z = - ω 2 ε 0 c 2 PNL,m , (18)
from where it is trivial to see that vector wave equation ( 16) is satisfied exactly, if the forward and backward components satisfy the following first order equations:

∂ Fm (z, ω) ∂z = - i 2k m (ω) ω 2 ε 0 c 2 PNL,m (z, ω)e +ikm(ω)z ∂ Bm (z, ω) ∂z = + i 2k m (ω) ω 2 ε 0 c 2 PNL,m (z, ω)e -ikm(ω)z . ( 19
)
It is worth noting that up to this point we did not make any assumptions, so the model is equivalent to the starting equations. Equations ( 19) represent a nonlinear boundary value problem that cannot be solved with direct methods, but must be solved iteratively. However in the great majority of cases of interest, we can assume that (i) there are no reflections and

(ii) that nonlinear polarization does not couple forward and backward waves (perturbative regime). In this case we can assume Bm (z, ω) ≈ 0 and Eqs. [START_REF] Kolesik | Unidirectional optical pulse propagation equation[END_REF], through [START_REF] Conforti | Modeling of ultrabroadband and single-cycle phenomena in anisotropic quadratic crystals[END_REF], reduce to the Forward Maxwell Equations:

∂ Êm (z, ω) ∂z + ik m (ω) Êm (z, ω) = -i ω 2ε 0 cn m (ω) PNL,m (z, ω). (20)
We consider an instantaneous nonlinear polarization composed of a quadratic and cubic parts (summation over repeated indexes is assumed)

P N L,j = ε 0 (χ (2) jkl E k E l + χ (3) jklm E k E l E m ), (21)
where χ

(2)

jkl and χ

jklm are the second and third order nonlinear susceptivity tensors, that are usually given in the crystal axes reference frame. In order to obtain the effective nonlinearity [START_REF] Midwinter | The effects of phase matching method and of uniaxial crystal symmetry on the polar distribution of second-order non-linear optical polarization[END_REF][START_REF] Midwinter | The effects of phase matching method and of crystal symmetry on the polar dependence of third-order non-linear optical polarization[END_REF], we have to rotate the polarization vector with matrix A, following the prescription

P N L (E) = AP N L (A T E). ( 22
)
After some calculations, we can write:

∂ Êx ∂z + ik x (ω) Êx = -iω cn x (ω) Px (23) ∂ Êy ∂z + ik y (ω) Êy = -iω cn y (ω)

Py

where the nonlinear terms P x , P y read as follows:

P x = d 0 E 2 x + 2d 1 E x E y + d 2 E 2 y + c 0 E 3 x + 3c 1 E 2 x E y + 3c 2 E 2 y E x + c 3 E 3 y , (24)
P y = d 1 E 2 x + 2d 2 E x E y + d 3 E 2 y + c 1 E 3 x + 3c 2 E 2 x E y + 3c 3 E 2 y E x + c 4 E 3 y . ( 25
)
where d m , m = 0, . . . , 3, are the effective nonlinearity for quadratic interactions, whereas c m , m = 0, . . . , 4 are the effective cubic nonlinearities. The values of the effective nonlinearity depend upon the crystal and their values can be found in literature [START_REF] Nikogosyan | Nonlinear Optical Crystals: A Complete Survey[END_REF][START_REF] Midwinter | The effects of phase matching method and of uniaxial crystal symmetry on the polar distribution of second-order non-linear optical polarization[END_REF][START_REF] Midwinter | The effects of phase matching method and of crystal symmetry on the polar dependence of third-order non-linear optical polarization[END_REF][START_REF] Banks | High intensity third-harmonic generation[END_REF]. In Tables I, II we report the effective nonlinearity for the crystals of class 3m, to which BBO belongs, and specify the kind of interaction. For example, eeo (e+e → o) indicates the sum frequency generation of two extraordinarily polarized electric fields (E x ) that generate an ordinarily polarized field (E y ).

Equations ( 23) are first order in the propagation coordinate, conserve the total field energy and retain their validity for arbitrary wide pulse bandwidth. The computational effort needed to solve these equations, by a standard split step Fourier method exploiting Runge-Kutta for the nonlinear step, is of the order of magnitude of that needed for solving the standard three-wave equations universally exploited to describe light propagation in quadratic crystals [START_REF] Baronio | Three-wave trapponic solitons for tunable high-repetition rate pulse train generation[END_REF][START_REF] Baronio | Velocity-locked solitary waves in quadratic media[END_REF]. However Eqs. ( 23) are far more general, and are equivalent to Maxwell equations when dealing with unidirectional propagation [START_REF] Kolesik | Unidirectional optical pulse propagation equation[END_REF][START_REF] Kolesik | Theory and simulation of ultrafast intense pulse propagation in extended media[END_REF].

III. EXAMPLES

In this section, we first show a representative example of the modeling of SHG of highintensity femtosecond pulses under strong influence of cubic nonlinearities. Then, we report quadratic ultrabroadband continuum dynamics with competing cubic nonlinearities. At last, we present soliton compression and dispersive waves dynamics in BBO, dominated by cubic nonlinearities.

A. High-intensity SHG

We fix the orientation angles of the BBO crystal to θ = 38 o and φ = 90 o .

We consider the propagation of an ordinary polarized (o) pulse of duration T = 30 fs, with intensity up to few TW/cm 2 , and central wavelength λ 0 = 630 nm injected in a 1mm long crystal. Under such assumptions, a type I (o + o → e) efficient phase-matched SHG interaction occurs.

Figure 1 shows the dependence of second harmonic conversion efficiency on the pump peak pulse intensity, with and without consideration of cubic nonlinearities. The obtained data indicate that for pump intensities exceeding 50 GW/cm 2 , the effects caused by cubic nonlinearities become significant and lead to a decrease of conversion efficiency. Figure 3a shows the time domain evolution of the ordinarily polarized (o) electric field envelope at 1200 nm during the propagation in BBO crystal. With envelope we mean the inverse Fourier transform of the positive frequency components of the spectrum. This visualization permit to have an envelope-like appearance, without fast oscillations of the carrier, but accounts of all frequency components. The input pulse undergoes a strong compression up to z = 0.6 mm, where the minimum pulse duration and maximum of spectral extension is achieved. Figure 3b shows the evolution of the ordinarily polarized field spectrum. The compression is due to the cascaded quadratic effects (γ 2 = -14 • 10 -16 m/V 2 , γ 3 = 6.65 • 10 -16 m/V 2 ). At the compression point the ordinary polarized pulse shows trailing oscillations, and subsequently radiation is emitted at a slower Figure 5a shows the time domain propagation of the ordinarily polarized electric field during the propagation in BBO crystal, figure 5b shows the evolution of the field spectrum.

The input pulse undergoes a strong compression up to z = 0.8 mm, where the minimum pulse duration and maximum of spectral extension. The compression is due to high-order cubic soliton excitation. A linear dispersive wave [START_REF] Wai | Soliton at the zero-group-dispersion wavelength of a single-model fiber[END_REF], located in the blue part of the spectrum at 900 nm, has been generated.

IV. CONCLUSIONS

We have derived unidirectional pulse propagation equations to describe extreme highintensity and ultra-broadband optical interactions in anisotropic crystals showing both quadratic and cubic nonlinear optical susceptibilities, taking BBO as the most relevant example. This model can be used to to enlighten high-order harmonic and ultrabroadband generation in BBO crystals under strong and competing influence of quadratic and cubic susceptivities.

FIG. 1 :

 1 FIG. 1: (Color online) Dependence of SHG efficiency on the pump peak intensity at the fundamental frequency, obtained without (blue stars) and with (red circles) consideration of the cubic nonlinear effects.

Figure 2

 2 Figure 2 shows typical evolutions of the field spectrum during the propagation in BBO crystal, at different intensity regimes. Above 50 GW/cm 2 , cubic nonlinearities give rise to a nonlinear-phase mismatch and to self-and cross-action of the interacting fundamental pump and second harmonics (i.e., o + o + o → o, o + o + e → e, e + e + o → o, e + e + e → e interactions), which not only decrease conversion efficiency but lead to spectral broadening of fundamental and second harmonic spectra.

FIG. 2 : 3 8ω

 23 FIG. 2: (Color online) Evolution of the total power spectrum | Êx | 2 + | Êy | 2 (decibels). The initial pulse has duration T = 30 fs, wavelength λ 0 = 630nm. The peak intensity is 10 GW/cm 2 (a), and 2 TW/cm 2 (b). Crystal's orientation: θ = 38 o and φ = 90 o .

  group velocity: a linear dispersive wave located in the red part of the spectrum at 2400 nm [? ].

FIG. 3 :

 3 FIG. 3: (Color online) a) temporal propagation and b) field spectrum evolution (decibels) of the ordinarily polarized electric field envelope in BBO crystal. The initial pulse has duration T = 20 fs, wavelength λ 0 = 1200 nm, peak intensity of 120 GW/cm 2 . Crystal's orientation: θ = 19 o and φ = 90 o .

Figure

  Figure4ashows the time domain evolution of the ordinarily polarized electric field during the propagation in BBO crystal, whereas figure4bshows the evolution of the field spectrum.The scenario has been dramatically changed with respect to the previous case. In fact, the effective quadratic negative Kerr nonlinearity (γ 2 = -6.6 • 10 -16 m/V 2 ), induced by mismatched type I (o + o → e) interaction, is perfectly balanced by the cubic nonlinearity of the medium (o + o + o → o interaction). The ordinarily polarized pulse propagates in the BBO crystal in the same way as the nonlinearities were vanishing, independently from input intensity.

FIG. 4 :

 4 FIG. 4: (Color online) a) temporal propagation and b) field spectrum evolution (decibels) of the ordinarily polarized electric field envelope in BBO crystal. The initial pulse has duration T = 20 fs, wavelength λ 0 = 1200 nm, peak intensity of 120 GW/cm 2 . Crystal's orientation θ = 16.2 o and φ = 90 o .

TABLE II :

 II Effective cubic nonlinear coefficients. c 11 = 5.6•10 -22 m 2 /V 2 , c 10 = -0.24•10 -22 m 2 /V 2 ,

	Coefficient		Expression	Interaction
	d 0	-3d 31 cos 2 θ sin θ -d 22 cos 3 θ sin 3φ	eee
	d 1		-d 22 cos 3φ cos 2 θ	eeo, oeo, oee
	d 2		-d 31 sin θ + d 22 cos θ sin 3φ	ooe, eoo, oeo
	d 3		d 22 cos 3φ	ooo
	TABLE I: Effective quadratic nonlinear coefficients. d 22 = 2.2pm/V, d 31 = 0.04pm/V [31].
	Coefficient	Expression	Interaction
	c 0		c 11 cos 4 θ + c 33 sin 4 θ + 3 2 c 16 sin 2 2θ	eeee
			-4c 10 sin 3φ sin θ cos 3 θ
	c 1		3 2 c 10 cos 3φ sin 2θ cos θ	eeoe, eeeo
	c 2		-1 3 c 11 cos 2 θ + c 16 sin 2 θ	ooee,eeoo
			+c 10 sin 2θ sin 3φ
	c 3		c 10 cos 3φ sin θ	oooe, ooeo
	c 4		c 11	oooo

c 16 = -1.4 • 10 -22 m 2 /V 2 [32].

V. ACKNOWLEDGMENTS We thank Daniele Faccio, Niclas Westerberg, and Stefano Trillo for discussions. The present research is supported in Brescia by the Italian Ministry of University and Research (MIUR), Grant PRIN 2009P3K72Z.


E Medernach
email: medernach@clermont.in2p3.fr

E Sanlaville
email: eric.sanlaville@math.univ-bpclermont.fr

Critère d'équité en ordonnancement de grille

Mesures pour l'équité Le problème de l'équité est un problème bien connu dans les réseaux telecoms. Les algorithmes de contrôle de congestion cherchent à allouer à chaque utilisateur un ux sur un graphe le plus équitablement possible en évitant un possible engorgement du réseau. Par exemple la norme ATM propose pour l'équité de maximiser le minimum des débits [START_REF] Bertsekas | Data Networks[END_REF]. Cependant le modèle du réseau est continu, et les ux calculés sont constants.

Un ordonnancement donné est caractérisé par son vecteur des vitesses d'exécution. Nous cherchons donc une mesure sur les vecteurs de R n + . Cette mesure peut être dénie comme une fonction de R n + dans R. Nous proposerons d'abord une dénition de l'équité via un certain nombre de caractéristiques désirables.

Denition 3 (Fonction équitable). Une fonction f symétrique est dite équitable si elle vérie :

∀i, j ≤ n, min(x i , x j ) < min(x i , x j ) ⇒ f (.., x i , .., x j , ..) ≤ f (.., x i , .., x j , ..) Les points de suspension indiquent des valeurs identiques entre elles. Si on cherche un ordre dans R n + au lieu d'une fonction scalaire on peut réecrire la dénition de l'équité, on trouve alors que le seul ordre qui vérie cette dénition est l'ordre lexicographique sur les vecteurs triés. Remarquons que cet ordre est un sur-ordre de l'ordre de majoration [START_REF] Bhargava | Using approximate majorization to characterize protocol fairness[END_REF] utilisé pour l'étude de protocoles réseaux.

Références

L

  'exigence de symétrie pour la fonction d'équité entraîne l'impartialité, cela exprime qu'aucun utilisateur n'est privilégié car f est invariante pour toute permutation des individus. Les indices peuvent être égaux ce qui signie que si on augmente l'allocation d'un seul utilisateur on augmente la fonction d'équité. Le critère d'équité est équivalent au fait que diminuer l'allocation à un utilisateur pour augmenter celle d'un autre déjà plus favorisé diminue l'équité. En eet si on passe de x à x soit la situation est dégradée pour chacun des utilisateurs soit la situation de celui qui a le moins empire alors que celle de l'autre s'améliore. On montre alors que toute fonction f continue qui vérie ces conditions est telle qu'il existe une fonction continue croissante g telle que f (x 1 , x 2 , . . .) = g(min(x 1 , x 2 , . . .)). En d'autres termes, une fonction vériant ces trois caractéristiques ne dépend que du minimum du vecteur d'entrée. Ceci conduit à un ordre partiel qui ne dépend que de la plus faible vitesse d'exécution. Cette conclusion n'est pas très satisfaisante.

1

  Contexte : la grille EGEE et la gestion des sites L'intérêt d'une grille de calcul est de permettre le partage de ressources de calcul et/ou de stockage. Dans un tel environnement concurrentiel apparaît le problème de distribuer ces ressources à la fois de façon ecace et équitable. Le contexte de notre travail est la grille EGEE qui se compose de 200 sites sur 50 pays pour un total de plus de 30000 CPUs. Un site d'EGEE se trouve au LPC-IN2P3 de Clermont Ferrand. L'infrastructure de grille EGEE est destinée aux applications 'on minimise le makespan on satisfait le système car l'ordonnancement obtenu est alors le plus compact possible. Mais on ne satisfait pas forcement les utilisateurs car certains peuvent être servis moins bien que d'autres. Ainsi minimiser une fonction économique globale qui ne prend pas en compte la présence de multiples utilisateurs peut aboutir à une situation injuste pour certains d'entre eux. Par exemple une politique favorisant les longs jobs comme LEPT (Longest Estimated Processing Time) donne de bons résultats pour le makespan mais pénalise les petits jobs, donc un utilisateur chez qui ce type de jobs prédomine. Critère individuel Chaque utilisateur utilise la puissance de calcul du site an d'accélérer ses temps de traitement. Il évalue donc la performance que lui propose un site par la vitesse de calcul accordée. L'objectif du site est de servir chaque utilisateur de façon équitable en orant une puissance de calcul bien répartie. Nous proposons de mesurer le niveau de service reçu individuel par la vitesse d'exécution moyenne obtenue de la part du site. Denition 2 (Vitesse d'execution). La vitesse d'exécution obtenue par l'utilisateur k est le ratio entre la quantité de calcul fournie au site par l'utilisateur k et la durée de traitement de cette quantité de calcul.

	Malheureusement à
	grande échelle ce choix crée un goulet d'étranglement et des problèmes de tolérance aux pannes.
	Au niveau global, il semble que la gestion entièrement centralisée des jobs entraîne de nombreux
	dysfonctionnements. En eet, l'absence de coordination entre sites ainsi qu'une latence importante
	des informations conduisent à de mauvais choix de placement des tâches.
	Au niveau local (gestion d'un site), il est fondamental pour une bonne gestion d'orir à chaque
	utilisateur en concurrence un accès équitable aux ressources. Par exemple dans les problèmes de les

de recherche scientique : les applications pilotes d'EGEE sont la physique des particules et les applications biomédicales. Plus d'une centaine de scientiques à travers le monde ont soumis environ 2 millions de jobs durant les dix premiers mois de l'année 2005. Ces jobs varient, de scripts très courts à des calculs s'étendant sur plusieurs jours.

Chaque utilisateur appartient à un groupe et est identié et connu par chaque site. L'objectif de la grille est de partager la puissance de calcul entre les utilisateurs ; le fait que les utilisateurs utilisent plusieurs sites permet de niveler la charge entre les sites. Les intergiciels de grille actuels ont recours à une architecture centralisée au dessus des moyens existants. d'attentes des études montrent qu'il est plus important d'être servi équitablement que de diminuer l'attente

[START_REF] Larson | Beyond the physics of queueing[END_REF]

.

Nous nous plaçons au cours de cette présentation au niveau local sur un site. Nous étudions comment ordonnancer les requêtes des utilisateurs de telle sorte que la part reçue par chacun ainsi que la part donnée aux groupes soit la plus équitable possible, dans un sens que nous cherchons à dénir.

2

Formalisation du problème étudié La soumission des tâches par utilisateur est soit irrégulière soit périodique pour certaines tâches spéciques. Le modèle d'ordonnancement sous-jacent est celui de tâches indépendantes exécutées sur des machines identiques en parallèle. L'ordonnancement se fait en ligne, le site reçoit les tâches au cours du temps sans connaître précisément leurs durées d'exécution. Celles-ci devraient être modélisées par des variables aléatoires diérentes suivant le groupe de l'utilisateur. L'utilisateur peut spécier une durée maximale pour ses tâches mais en pratique cette mesure se révèle imprécise. Denition 1. Le problème mono-site et mono-groupe est donc déni par : m machines identiques n utilisateurs p j,k : durée du j e job de l'utilisateur k, r j,k : date d'arrivée. Les jobs sont supposés monoprocesseurs et indépendants les uns des autres. L'objectif est de maximiser l'équité. Selon des études préliminaires [2], le temps entre deux arrivées sur les clusters locaux du LPC peut être modélisé par une distribution hyper-exponentielle par utilisateur. Les durées des tâches varient beaucoup (de quelques secondes à plusieurs jours) et suivent une loi log-linéaire par morceaux. La notion d'équité LorsquLa vitesse d'exécution peut aussi être appelée débit de calcul. Clairement, la diérence face aux critères classiques est la prise en compte de l'utilisateur considéré comme un sous-ensemble de jobs indépendants.