V Denis

A Pelat

C Touzé

F Gautier

Improvement of the acoustic black hole eect by using energy transfer due to geometric nonlinearity

Acoustic Black Hole eect (ABH) is a passive vibration damping technique without added mass based on exural waves properties in thin structures with variable thickness. A common implementation is a plate edge where the thickness is locally reduced with a power law prole and covered with a viscoelastic layer. The plate displacement in the small thickness region is large and easily exceeds the plate thickness. This is the origin of geometric nonlinearity which can generate couplings between linear eigenmodes of the structure and induce energy transfer between low and high frequency regimes. This phenomenon may be used to increase the eciency of the ABH treatment in the low frequency regime where it is usually inecient. An experimental investigation evidenced that usual ABH implementation gives rise to measurable geometric nonlinearity and typical nonlinear phenomena. In particular, strongly nonlinear regime and wave turbulence are reported. The nonlinear ABH beam is then modeled as a von Kármán plate with variable thickness. The model is solved numerically by using a modal method combined with an energy-conserving time integration scheme. The eects of both the thickness prole and the damping layer are then investigated in order to improve the damping properties of an ABH beam. It is found that a compromise between the two eects can lead to an important gain of eciency in the low frequency range.

Introduction

Controlling vibration is a major concern in many industrial applications today. Classical methods for passive vibration damping include the use of heavy viscoelastic layer [START_REF] Ross | Damping of plate exural vibrations by means of viscoelastic laminae[END_REF] or tuned mass dampers [START_REF] Jones | The eect of tuned dampers on vibrations of simple structures[END_REF]. These techniques are ecient and largely studied but their implementation is often limited by an important added mass or a narrow frequency range. The Acoustic Black Hole (ABH) eect is a particular surface damping method adapted to mitigate vibrations on a wide frequency range without added mass. The name ABH in this context has been rst introduced by Mironov [START_REF] Mironov | One-dimensional acoustical "black hole[END_REF] and Krylov [START_REF] Krylov | Acoustic 'black holes' for exural waves as eective vibration dampers[END_REF] and is now used in an increasing number of references in the literature [START_REF] Georgiev | Damping of structural vibrations in beams and elliptical plates using the acoustic black hole eect[END_REF][START_REF] Zhao | Experimental verication of energy harvesting performance in plate-like structures with embedded acoustic black holes[END_REF][START_REF] Huang | Flexural wave focalization in plates with imperfect two-dimensional acoustic black hole[END_REF]. Although the name may not be completely adapted as referring to an astrophysical phenomenon it is kept in this paper in order to be consistent with the current developments around this particular vibration absorber.

As described by Mironov [START_REF] Mironov | Propagation of a exural wave in a plate whose thickness decreases smoothly to zero in a nite interval[END_REF], the ABH eect takes advantage of the exural wave properties in a plate with variable thickness. It can be shown that waves propagating in a beam or plate with thickness h(x)

h(x) = x m , m ≥ 2, (1)
slow down without being reected and that their travel time in the prole tends to innity if thickness tends to zero. Due to the small thicknesses involved, manufacturing this "Acoustic Black Hole" prole is dicult and the residual thickness at the edge may be too large to obtain the desired eect. The addition of a thin viscoelastic layer [START_REF] Krylov | Acoustic 'black holes' for exural waves as eective vibration dampers[END_REF] can compensate the eect of the prole truncation and the reection coecient of the extremity is then considerably reduced. Further theoretical works propose a more rened model for the reection coecient [START_REF] Georgiev | Damping of structural vibrations in beams and elliptical plates using the acoustic black hole eect[END_REF] of an ABH termination and demonstrate the resulting increase of the modal overlap [START_REF] Denis | Modal overlap factor of a beam with an acoustic black hole termination[END_REF]. Experimental evidences of the eect are numerous [START_REF] O'boy | Damping of exural vibrations in rectangular plates using the acoustic black hole eect[END_REF][START_REF] Bowyer | Eect of geometrical and material imperfections on damping exural vibrations in plates with attached wedges of power law prole[END_REF][START_REF] Denis | Measurement and modelling of the reection coecient of an acoustic black hole termination[END_REF]. While it is very ecient at high frequencies, several works point out the ineciency of the ABH technique in the low frequency regime [START_REF] Denis | Modal overlap factor of a beam with an acoustic black hole termination[END_REF][START_REF] Conlon | Numerical analysis of the vibroacoustic properties of plates with embedded grids of acoustic black holes[END_REF], which may be a drawback in its potential applications.

The minimum thickness of the ABH extremity needed to obtain an important damping is very small (about hundred microns) and the displacement amplitudes at the extremity become large with respect to the thickness. Structures presenting this behavior are usually studied in the framework of nonlinear vibrations [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. To the authors' knowledge, all numerical and theoretical models of ABH are based on linear assumptions such that the eects of the large amplitudes are neglected. Nonlinear eects include couplings between linear eigenmodes of the structure [START_REF] Nayfeh | Nonlinear oscillations[END_REF][START_REF] Touzé | Transition scenario to turbulence in thin vibrating plates[END_REF]. Such couplings could be used for transferring energy from low frequency regime, where the ABH is inactive, to high frequency regime where vibrations can be eciently dissipated, thus improving the ABH performance. This paper presents an investigation of the nonlinear behavior of ABH beams and evaluates the interest of these nonlinear eects to improve the ABH damping eciency at low frequency. Sec. 2 reports an experimental highlight of nonlinear response of an ABH beam. Sec. 3 develops a von Kármán plate model with variable thickness for the ABH beam and its numerical resolution from a modal projection and an energy-conserving time-marching scheme. The results of the model are detailed in Sec. [START_REF] Krylov | Acoustic 'black holes' for exural waves as eective vibration dampers[END_REF]. Conclusions and perspectives are given in Sec. 5.

Experiment

This section is devoted to an experimental investigation of the nonlinear behavior of a typical ABH beam. Dedicated forced vibration experiments are conducted in order to reveal the nonlinearities which can induce a variety of typical nonlinear phenomena including wave turbulence, characterized by a broadband Fourier spectrum and an energy transfer from low to high frequency range.

Experimental setup

The experimental study concerns two aluminium beams : the rst is uniform with dimensions 1.3 m × 20 mm × 5 mm. Given its thickness, this beam is awaited to behave linearly. The second selected beam, with dimensions 1.5 m × 20 mm × 5 mm, has an Acoustic Black Hole termination of length 0.35 m; the thickness of this extremity decreases quadratically down to 50 µm approximately (see the geometry of the sample in Fig. 3).

A thin viscoelastic tape (for example an electric tape) is applied on the ABH extremity. It is approximately 50 µm thick and its length is 0.2 m.

The experimental setup consists in a shaker (LDS V201), an amplier (LDS PA25E), a force sensor (PCB 208C03), an accelerometer (PCB 352C29) and an acquisition card (NI-USB4431). It is driven by a Labview interface and the results are exported and processed with MATLAB. The beam to be tested is vertically suspended and excited at its midpoint with the shaker. The accelerometer is placed at an arbitrary location on the beam in order to record the vibrations at a given point.

Nonlinear regimes

The reference and ABH beams are excited at 83 Hz with a forcing amplitude increasing from 0 to 17 N in 40 s. By increasing gradually the forcing amplitude, the transition from periodic to chaotic regime can be observed [START_REF] Touzé | Transition to chaotic vibrations for harmonically forced perfect and imperfect circular plates[END_REF], and compared to the transition scenario given in [START_REF] Touzé | Transition scenario to turbulence in thin vibrating plates[END_REF]. Figs. 1(a-b) show spectrograms of output acceleration for the two congurations. The spectrum of the response for the uniform beam shows that most of the energy is concentrated in the excitation frequency at 83 Hz, which is an indication that the beam behaves linearly (see Fig. 1(a)). Additional components in the spectrum may be observed. Their amplitudes are small compared to the peak at the excitation frequency. They can be attributed to other elements in the experimental setup, or residual nonlinearities in the shaker. The behavior is signicantly dierent for the ABH beam, see Fig. 1(b), the spectrum of which rapidly enriches with a very large number of harmonics present in the response. The response of the beam is periodic but nonlinear, with a spread of energy in a large harmonic series which, in this reported case, goes up to 3 kHz. Similar results have been obtained with numerous dierent forcing frequencies. Assuming that the beam's material constitutive behavior is linearly elastic, the nonlinearity present in the system comes from large amplitude oscillations of the beam, with respect to the thickness, and is denoted as geometric nonlinearity in the literature [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. It signicantly inuences the system's response. Fig. 1(c) has been obtained for a dierent forcing frequency, namely 113 Hz. It reveals the case of a complete transition to the turbulent, chaotic behavior. A direct transition from the periodic regime, composed of numerous harmonics, to the turbulent regime, characterized by a broadband Fourier spectrum and a cascade of energy to the high frequencies [START_REF] During | Weak turbulence for a vibrating plate: can one hear a Kolmogorov spectrum?[END_REF][START_REF] Boudaoud | Observation of wave turbulence in vibrating plates[END_REF][START_REF] Cadot | Wave turbulence in vibrating plates[END_REF], is observed after 30 seconds of increasing the forcing amplitude. This observation is in the line of the transition to turbulence scenario for thin vibrating plates, already commented and analyzed in details with experiments and numerical simulations in [START_REF] Touzé | Transition scenario to turbulence in thin vibrating plates[END_REF][START_REF] Touzé | Transition to chaotic vibrations for harmonically forced perfect and imperfect circular plates[END_REF]. A more peculiar feature is shown in Fig. 1(d), where the ABH beam excited at 102 Hz rst transits to the turbulent regime, clearly observed between 10 and 20 seconds, and then, while the forcing amplitude still increasing, recovers a periodic motion, characterized by a Fourier spectrum composed of discrete peaks. Note that the regime is still highly nonlinear as the number of harmonics is very large with energy up to 8 kHz. This peculiar regime has been observed a signicant number of times during all the experiments realized, so that it appeared meaningful to report this behavior even if it is not analyzed specically. During these large amplitude periodic nonlinear oscillations, it has been observed that the ABH beam vibrates around a deected position, dierent from the position of the beam at rest. Hence these oscillations can be interpreted as strongly nonlinear normal modes (NNM), i.e. periodic orbits, existing at large energy levels, and far from the rest (straight) position. They can be brought closer to earlier experimental observations shown on vibrating beams in [START_REF] Cusumano | Chaotic non-planar vibrations of the thin elastica, part I: experimental observation of planar instability[END_REF]. In these strongly nonlinear regimes, the system may have a number of large-amplitude NNMs that are isolas in the phase space, and for some particular conditions, nonlinear oscillations can be trapped in one of these regimes, as observed in Fig. 1(d).

Energy transfer between low and high frequencies

In order to better understand the nonlinearity and the benet that can be gained from it, the beams are now excited with a ltered white noise in the 100-500 Hz range, and for several increasing amplitudes. This frequency band has been selected as being the band in which the ABH is generally not ecient enough, and where the nonlinearity can bring an improvement by transferring the energy to higher frequencies. The parameter σ denotes the standard deviation of the exciting force for each case. For a white noise, σ is proportional to the square root of the power. The Power Spectral Density (PSD) of the acceleration a(t), normalized by σ 2 , is dened as:

S aa (f ) = 1 σ 2 ∞ -∞ ∞ -∞ a * (t)a(t -τ )e -j2πf τ dt dτ. (2)
It is displayed in Fig. 2 for uniform and ABH beams. Fig. 2(a) shows that the uniform beam exhibits a linear behavior since the power is concentrated in the excitation range, except for few artefacts that may come from the setup and happen however at a much lower scale. Most importantly, the spectrum does not depend on the amplitude. On the contrary, in the ABH case (see Fig. 2(b)), the power largely leaks outside the excitation range for high levels of excitation. On the zoom of Fig. 2(c) a small reduction of the vibratory levels with the increase of amplitude can actually be observed for the ABH beam, more specically on the two peaks around 360 and 450 Hz. This is interpreted as an energy transfer from low to high frequency domains, the latter being inherently more damped. This is particularly interesting as it could be used in order to increase the eciency of the ABH treatment at low frequencies. Note that here the transfer is not important enough and the gain in the damping of vibrations is very limited. The aim of the present study is thus to develop a model of the ABH including geometric nonlinearity, in order to analyze how this transfer can be optimized for a better eciency of the ABH in the low frequency range.

3 Model of nonlinear acoustic black hole

A von Kármán plate model with variable thickness

The ABH beam is modelled as a nonlinear von Kármán plate [START_REF] Nayfeh | Nonlinear oscillations[END_REF][START_REF] Von Kármán | Festigkeitsprobleme im maschinenbau[END_REF] with variable thickness. It is emphasized that the plate prole is assumed to vary along the x direction only in this study, hence plate characteristics such as mass density ρ and thickness h only depend on x. The equation of motion for the transverse displacement w(x, y, t) can then be written, following the developments proposed in [START_REF] Efstathiades | A new approach to the large-deection vibrations of imperfect circular disks using Galerkin's procedure[END_REF][START_REF] Thomas | Geometrically nonlinear exural vibrations of plates: In-plane boundary conditions and some symmetry properties[END_REF], as:

ρ(x)h(x) ẅ + ∆(D(x)∆w) -(1 -ν)L(D(x), w) = p + L(w, F ), (3)
where D(x), ρ(x), h(x) and ν are respectively bending stiness, mass density, thickness and Poisson ratio of the plate. The term p(x, y, t) represents the exciting force per unit surface. The Airy stress function F (x, y, t) respects the following compatibility equation:

∆( 1 E(x)h(x) ∆F ) -(1 + ν)L( 1 E(x)h(x) , F ) = - 1 2 L(w, w), (4)
where E(x) is the Young modulus of the plate. Finally L(f, g) is the Monge-Ampère bilinear operator which can be written, in cartesian coordinates, as:

L(f, g) = f xx g yy + f yy g xx -2f xy g xy , (5)
where f xx denotes the second derivative of f with respect to x.

The geometry of the plate is shown in Fig. 3. It is composed of three dierent regions. The rst one on the right is a uniform region (primary structure) with constant thickness h 0 . The second region is the ABH where the thickness decreases from h 0 to h t , which represents the minimum thickness of the plate. An extension of the small thickness region is considered in the left part of the beam, which is parameterized by its length l add that can be tuned in order to increase the nonlinear behavior of the structure. The thickness h p (x) of the plate can thus be written as: where m is the exponent of the ABH prole and l ABH is the length of the region of variable thickness. The plate has a width b and a total length L. In all the cases tested in this study, the exponent m is selected as m=2.

h p (x) =        h t ∀x ≤ l add , h t + (h 0 -h t ) (x-l add ) m l m ABH ∀x ∈ [l add , l add + l ABH ], h 0 ∀x ≥ l add + l ABH , (6)
The eect of the damping layer is represented with equivalent mechanical properties that modify the bending stiness, the mass density and the thickness, following the model of Ross-Ungar-Kerwin [START_REF] Ross | Damping of plate exural vibrations by means of viscoelastic laminae[END_REF]. As shown in Fig. 3, the added viscoelastic layer is assumed to be of thickness h l , and located in the interval [x l1 , x l2 ]. Its eect on the exural rigidity is taken into account more easily in the frequency domain, where the damping is represented as an imaginary part of the complex bending stiness D * (x), which can be written as:

D * (x) =                    E p h p (x) 3 12(1 -ν 2 ) (1 + jη p ) = D p (x)(1 + jη p ), ∀x / ∈ [x l1 , x l2 ], D p (x) (1 + jη p ) + E l Ep h l h p (x) 3 (1 + jη l ) + 3 1 + h l hp(x) 2 E l h l Ephp(x) (1 -η p η l + j(η p + η l )) 1 + E l h l Ephp(x) (1 + jη l ) , ∀x ∈ [x l1 , x l2 ], (7)
where j is the imaginary unit, D p (x), E p , η p and h p (x) are respectively the bending stiness, the Young modulus, the loss factor and the thickness of the plate without the layer, while E l , η l and h l are the Young modulus, the loss factor and the thickness of the damping layer, respectively. The real bending stiness D(x) is dened as the real part of D * (x):

D(x) = Re (D * (x)) , (8)
which is nally inserted into the equations of motion ( 3)-( 4) in order to take into account all these eects in the nal model. The eect of the mass added by the presence of the viscoelastic layer is taken into account by introducing the following equivalent mass density ρ(x):

ρ(x) =    ρ p , ∀x / ∈ [x l1 , x l2 ], ρ p h p (x) + ρ l h l h p (x) + h l , ∀x ∈ [x l1 , x l2 ], (9)
where ρ p and ρ l are the mass densities of the structure and the layer respectively. The total thickness h(x) of the sandwich plate can nally be written as:

h(x) = h p (x), ∀x / ∈ [x l1 , x l2 ], h p (x) + h l , ∀x ∈ [x l1 , x l2 ]. (10)
The equivalent Young modulus of the sandwich plate that have to be inserted in (4) for the compatibility equation on the Airy stress function, is deduced from the real bending stiness D(x) and the total thickness as:

E(x) = 12D(x)(1 -ν 2 ) h(x) 3 . (11)
x z

x l 2 x l 1 0 h 0 l ABH l add h t h l h p (x) = h t + (x -l add ) m (h 0 -h t )/l m ABH ⊗ y
Figure 3: Schematic view of the ABH plate thickness prole. h 0 is the thickness of the plate thickness in the uniform region, h t is the minimum thickness of the plate, m is the exponent of the ABH prole, l ABH is the length of the region of variable thickness and l add is the length of the extension of small constant thickness. The damping layer stands between abscissas x l1 and x l2 .

More concisely, the equations of motion can be written as:

ρh ẅ + (D, w) = p + L(w, F ), (12a)
♦(A, F ) = -L(w, w)/2, (12b)
with A(x) = 1/E(x)h(x), and and ♦ the spatial operators applied to w and F :

(D, w) = ∆(D(x)∆w) -(1 -ν)L(D(x), w), (13)
♦(D, w) = ∆( 1 E(x)h(x) ∆F ) -(1 + ν)L( 1 E(x)h(x) , F ). (14)
In the case of a plate with free edges, the boundary conditions write, on an edge with normal line n and tangent one τ [START_REF] Denis | Modal overlap factor of a beam with an acoustic black hole termination[END_REF][START_REF] Thomas | Geometrically nonlinear exural vibrations of plates: In-plane boundary conditions and some symmetry properties[END_REF]:

- ∂ ∂n D ∂ 2 w ∂n 2 + ν ∂ 2 w ∂τ 2 -2(1 -ν) ∂ ∂τ D ∂ 2 w ∂n∂τ = 0, (15a)
∂ 2 w ∂τ 2 + ν ∂ 2 w ∂n 2 = 0, (15b)
F = ∂F ∂n = 0. ( 15c
)
One can note that Eq. (15c) are not the classical free conditions but are sucient conditions to impose an edge free of loads in the in-plane direction [START_REF] Thomas | Geometrically nonlinear exural vibrations of plates: In-plane boundary conditions and some symmetry properties[END_REF].

Modal approach

A modal approach is used for the spatial discretization of the problem, following the general procedure described in [START_REF] Ducceschi | Modal approach for nonlinear vibrations of damped impact plates: application to sound synthesis of gongs and cymbals[END_REF]. Here the diculty comes from the plate's variable thickness which, in turn, gives rise to a more complex formulation for the linear operators given in Eqs. [START_REF] Denis | Measurement and modelling of the reection coecient of an acoustic black hole termination[END_REF]. The eigenmodes are computed for the unforced, conservative problem. For that purpose, the real bending stiness D, dened in ( 8) is used. The conservative and linear associated problems to be solved nally read:

-ρhω 2 w + (D, w) = 0, (16)
and

♦(A, F ) = ζ 4 F, (17)
to which are added the boundary conditions [START_REF] Touzé | Transition scenario to turbulence in thin vibrating plates[END_REF]. These are eigenvalue problems with solutions (ω k , Φ k (x, y)) and (ζ n , Ψ n (x, y)), respectively. The displacement and the Airy stress function are then expressed as:

w(x, t) = NΦ k=1 Φ k (x, y)q k (t), (18)
and

F (x, t) = NΨ n=1 Ψ n (x, y)β n (t), (19)
where q k (t) and β n (t) are the modal coordinates for the transverse displacement and the Airy stress function, respectively. The integers N Φ and N Ψ represent the number of transverse and in-plane modes kept in the modal truncation.

It is emphasized that the modes are orthogonal with respect to the mass:

S ρ(x)h(x)Φ k (x, y)Φ p (x, y)dS = 0, p = q, (20)
and normed such that:

S Φ k (x, y) 2 dS = 1. (21)
The projection on the linear modes gives two equivalent formulations of the problem. The formulation ( 22) is quadratic in (q, β) and can be written as:

β p (t) = - 1 2ζ 4 p NΦ i=1 NΦ j=1 H p i,j q i (t)q j (t), (22a)
qs (t) + ω 2 s q s (t) = F s (t) + 1 M s NΦ k=1 NΨ n=1 E s k,n q k (t)β n (t), (22b)
with M s = S ρ(x)h(x)Φ s (x, y) 2 dS the modal masses, F s = 1 Ms S pΦ s dS the modal force, and the coupling coecients:

H p i,j = S Ψ p L(Φ i , Φ j )dS, (23)
and

E s k,l = S L(Φ k , Ψ l )Φ s dS. (24)
The formulation [START_REF] Leissa | Vibration of plates[END_REF], where variables β have been eliminated to keep only the modal amplitudes q of the transverse displacement w, shows a cubic nonlinearity in q:

qs (t) + ω 2 s q s (t) = F s (t) - 1 M s NΦ k=1 NΦ m=1 NΦ n=1 Γ s k,m,n q k (t)q m (t)q n (t). (25)
The Γ s k,m,n coecients are here calculated from H l m,n and E s k,l :

Γ s k,m,n = NΨ l=1 H l m,n E s k,l 2ζ 4 l . ( 26
)
Although completely equivalent, the two dierent formulations can be used for dierent purposes. As shown in section 3.5, an energy-conserving scheme can be adopted for a robust and accurate time integration of the system. As shown in [START_REF] Ducceschi | Modal approach for nonlinear vibrations of damped impact plates: application to sound synthesis of gongs and cymbals[END_REF], this scheme particularly suits the quadratic formulation. On the other hand, the cubic formulation [START_REF] Leissa | Vibration of plates[END_REF] shows that in-plane variables β can be eliminated from the procedure. As two dierent integers have been introduced in the modal truncation given in Eqs. ( 18)-( 19), one can use the cubic formulation in order to verify the convergence of the solution as function of the number of in-plane modes N Ψ . Indeed, monitoring the convergence of the terms of the tensor Γ s k,m,n versus N Ψ gives the sucient number of in-plane modes needed for convergence. This will be shown in more details in the result section, see 4.2.

One can remark that the properties H p i,j = H p j,i and H p i,j = E i p,j already described in [START_REF] Thomas | Geometrically nonlinear exural vibrations of plates: In-plane boundary conditions and some symmetry properties[END_REF][START_REF] Ducceschi | Modal approach for nonlinear vibrations of damped impact plates: application to sound synthesis of gongs and cymbals[END_REF] for the coupling tensors are also valid in this case and allow one to reduce the number of computed coecients and thus the computational burden associated with this preliminary step of oine calculations.

Computation of the linear modes

Due to the variable thickness, the solution to the eigenvalue problems [START_REF] Touzé | Transition to chaotic vibrations for harmonically forced perfect and imperfect circular plates[END_REF][START_REF] During | Weak turbulence for a vibrating plate: can one hear a Kolmogorov spectrum?[END_REF] is particularly dicult and needs a dedicated numerical approach. The solutions are obtained using a nite dierence method thoroughly described in [START_REF] Denis | Modal overlap factor of a beam with an acoustic black hole termination[END_REF]. The use of a uniform grid for a structure where the wavelength varies spatially can result in a high computational cost. This problem is solved using a coordinate change and a transformed problem, which is equivalent to using a grid adapted to the exural wavelength variation. Since the wavelength λ = 2π 4 Eh 2 /12ρ(1ν 2 )ω 2 varies as the square root of the thickness according to classical plate theory [START_REF] Leissa | Vibration of plates[END_REF], the coordinate change is selected as follows:

dx = dx xav h(x) , with xav = L 0 1 h(x) dx, (27)
where x denes the transformed coordinate.

First the displacement w(x, y), the bending stiness D(x), the mass density ρ(x) and the thickness h(x) are approximated with grid functions w n,q , D n , ρ n and h n , respectively, where n and q are the spatial indexes such that x = n∆ x and y = q∆ y , with ∆ x and ∆ y the spatial steps along the x and y directions. Then the eigenvalue problem in transverse displacement ( 16) is approximated with a transformed discrete problem, which can be written as [START_REF] Bilbao | Numerical sound synthesis : nite dierence schemes and simulation in musical acoustics[END_REF][START_REF] Denis | Modal overlap factor of a beam with an acoustic black hole termination[END_REF]:

h -1/2 n x2 av δ x+ (µ x-h -1/2 n )δ x-D n h -1/2 n x2 av δ x+ ((µ x-h -1/2 n )δ x-w n,q ) + δ yy w n,q + D n δ yy h -1/2 n x2 av δ x+ ((µ x-h -1/2 n )δ x-w n,q ) + δ yy w n,q -(1 -ν) h -1/2 n x2 av δ x+ ((µ x-h -1/2 n )δ x-D n )δ yy w n,q -ρ n h n ω 2 w n,q = 0. ( 28
)
where δ x-, δ x+ , µ x-, µ x+ and δ yy are the nite dierence operators dened in A.

The in-plane problem is treated in the same manner. The convergence of the method has been extensively studied in [START_REF] Denis | Modal overlap factor of a beam with an acoustic black hole termination[END_REF]. The results presented in this article have been obtained using a 2000 × 100 grid. It is reminded that nite dierence methods (especially those using second-order approximations), although converged, still presents an inherent numerical dispersion [START_REF] Bilbao | Numerical sound synthesis : nite dierence schemes and simulation in musical acoustics[END_REF]. Other numerical methods, such as Finite Element Analysis (FEA), are indeed also suitable to nd the desired solutions, as it has for instance been recently done in the current context in [START_REF] Conlon | Numerical analysis of the vibroacoustic properties of plates with embedded grids of acoustic black holes[END_REF] for circular ABH.

The resolution of the discrete eigenvalue problems yields the out-of-plane and in-plane modes and the associated eigenvalues. The coecients H p i,j , E s k,l and Γ s k,m,n are computed with the expressions (23, 24, 26) from the numerical mode shapes, following the general procedure also described in [START_REF] Touzé | Direct nite element computation of non-linear modal coupling coecients for reduced-order shell models[END_REF] for nite element shell models.

Losses

In accordance with the general modal framework used to solve the problem, the losses are taken into account by using a modal damping ratio for each oscillator equation. Assuming they remain small, the modal damping ratios ξ s are obtained by solving the eigenvalue problem in displacement with complex bending stiness:

-ρhω * 2 w + (D * , w) = 0. ( 29
)
According to the literature [START_REF] Géradin | Mechanical Vibrations: Theory and Applications to Structural Dynamics[END_REF], the eigenvalues ω * k of the problem dened in Eq. ( 29) can be written as:

ω * s = -ω s ξ s -j 1 -ξ 2 s , (30)
where ωs are the natural frequencies and ξ s are the modal damping ratios. The modal damping ratios are then directly used in the system (22a-22b), which can nally be written as:

β p (t) = - 1 2ζ 4 p NΦ i=1 NΦ j=1 H p i,j q i (t)q j (t), (31)
qs (t) + 2ξ s ω s qs (t) + ω 2 s q s (t) = F s (t) + 1 M s NΦ k=1 NΨ n=1 E s k,n q k (t)β n (t). ( 32
)
This method is advantageous since it yields a dissipative time evolution problem that can be solved in the same manner as for the conservative one. Note that modal damping have been introduced only on Eq. ( 32) for the transverse motions w. This can be justied by considering the cubic formulation given in Eq. [START_REF] Leissa | Vibration of plates[END_REF]. Indeed the in-plane motions can be condensed in order to use a closed formulation for the modal displacements q s only where damping has to be introduced. Secondly, inertia eects are not taken into account for the in-plane motions, which gives a second justications for neglecting possible losses in Eq. (31).

Time integration scheme

The system of Eqs. [START_REF] Ostiguy | Eects of initial geometric imperfections on dynamic behaviour of rectangular plates[END_REF][START_REF] Ducceschi | Dynamics of the wave turbulence spectrum in vibrating plates: A numerical investigation using a conservative nite dierence scheme[END_REF] is solved using an energy-conserving scheme described in [START_REF] Bilbao | A family of conservative nite dierence schemes for the dynamical von Kármán plate equations[END_REF][START_REF] Ducceschi | Modal approach for nonlinear vibrations of damped impact plates: application to sound synthesis of gongs and cymbals[END_REF], that is conditionally stable. The rst step is to discretize the modal coecients q s (t) → q n s (where t = n∆t and ∆t is the time step) and the derivation operators with the proper nite dierence operators. The conservative scheme can then be written [START_REF] Ducceschi | Modal approach for nonlinear vibrations of damped impact plates: application to sound synthesis of gongs and cymbals[END_REF]:

δ tt q n s + 2ξ s ω s δ t. q n s + ω 2 s q n s = 1 M s NΦ k=1 NΨ l=1 E s k,l q n k [µ t. β n l ] + F n s , (33)
µ t-β n p = -1 2ζ 4 p NΦ i=1 NΦ j=1 H p i,j q n i [e t-q n j ], (34)
where the operators δ tt , δ t. , µ t-, µ t. and e t-are dened in A.

The stability condition of such a scheme is [START_REF] Ducceschi | Modal approach for nonlinear vibrations of damped impact plates: application to sound synthesis of gongs and cymbals[END_REF]:

f e > πf max , (35)
where f e = 1/∆t is the sampling rate and f max is the largest frequency retained in the modal truncation for the transverse motion w, i.e. f max = f NΦ .

4 Numerical results

Parameters of the simulated beams

The model is used to analyze and compare the damping performance of thirteen dierent beams, denoted U, E0, E0D0, E1, E2, E3, E4, E4D0, E4D1, E4D2, E4D3 and E4D4, and represented schematically in Fig. 4.

A special emphasis is put on their ability to improve the eciency of the ABH in the low frequency regime. All the tested beams have common parameters that are gathered in Tab. 1, while the distinctive features of each of these congurations are given in Tab. 2. Conguration U is a uniform beam and will be used as a reference case for comparing the dierent eects brought respectively by the ABH, the nonlinear additional length and the damping layer. Beams E0 to E4 are ABH beams without viscoelastic layer. Their additional length l add is gradually increased in order to extend the interval of the beam where the thickness is minimal. Their behavior is expected to be more and more nonlinear. Conguration E0D0 has the same geometry than E0, the only dierence being the presence of the viscoelastic layer. Finally congurations E4D0 to E4D4 are all based on E4, with the presence of the damping layer of increasing length. As a memory aid, E stands for extension, D stands for damping and the number immediately after the letter is an indication for the concerned length.

Table 1: Values of the parameters shared by all the congurations simulated with the numerical model.

Parameters of the beam

Ep [GPa] ρp [kg m -3 ] ν ηp [%] h 0 [mm] L [m] b [
E l [GPa] ρ l [kg m -3 ] η l [%] h l [µm]
7 1000 40 100

Convergence study

The convergence of the numerical scheme is rst investigated. The main free parameters that have to be set are the number of modes kept in the truncation, both for transverse and in-plane modes, N Φ and N Ψ , and the time step ∆t, or, equivalently, the sampling rate f e . In order to ensure a converged representation of the in-plane motions, the convergence of coecient Γ p p,p,p with respect to N Ψ is investigated. An example is shown in Fig. 6, for p=29 (which correspond to the transverse mode (23,0) with 23 nodal lines along x and no nodal lines along y, the modes being sorted by increasing eigenfrequencies), for two dierent congurations : the reference beam U and the ABH E0 with a small length of decreasing thickness and no damping layer. The corresponding modal shape for E0 is plotted in Fig. 5(a).

In all cases, the coupling term Γ p p,p,p for the simple conguration U converges very quickly: approximately 100 in-plane modes are sucient in order to ensure a converged value and thus a ne representation of the inplane motions. On the other hand, the convergence for the beam E0 is much slower, as evidenced in Fig. 6(a). Note also that the value of the coupling coecient is also much larger in case E0 than U, which was awaited since the small thickness region is designed in order to increase the nonlinearity in the system. The rst consequence is thus the values of the nonlinear coupling coecients, which are substantially increased. The convergence of Γ p p,p,p for mode p=29 is slow and behaves incrementally; which suggests that some in-plane modes contributes strongly to the nal value of the coecients while others do not. In fact, the convergence can be made much faster by sorting the in-plane modes by families instead of sorting them by increasing eigenfrequencies. Let us dene in-plane family n the family containing the modes with n+1 nodal lines across the y-dimension (including the boundaries). Two examples are shown in Fig. 5(c-d): mode Ψ 101,2 in Fig. 5(c) belongs to family 1 while mode Ψ 61,3 in Fig. 5(d) belongs to family 2. Fig. 6(b) shows the result of the computation when the in-plane modes are sorted by families. Remarkably enough, only modes of family 1 actually contribute to the value of the coupling term for this specic transverse mode p=29. The same remarks can be made for other value of p and it is veried that these result hold for N Φ = 200 out-of-plane modes in the modal truncation. The sorting procedure allows to eliminate non-contributing families from the computation, which greatly reduces the computational cost. The number of transverse modes N Φ depends on the frequency band one would like to cover in the dynamical simulations. It depends on the nonlinear regime at hand and the largest frequency present in the spectrum of the solution. For all the presented simulations, N Φ =200 has been selected, which gives an accurate representation of the solution up to the frequency 7 kHz for the most nonlinear conguration E4. This number has been found to be sucient for the tested regimes. Finally, once the number of in-plane modes chosen, the sampling rate f e can be simply selected thanks to the stability rule given by Eq. [START_REF] Miquel | Nonlinear dynamics of exural wave turbulence[END_REF]. However a larger value has been retained with f e =100 kHz in all the computations, in order to have a more rened time representation.

Nonlinear regimes

In this section, the behavior of the model is investigated with respect to a harmonic forcing:

p(x, y, t) = δ(x -x 0 )δ(y -y 0 )G(t) sin(Ωt), (36)
localized at point (x 0 ,y 0 ), with xed frequency Ω and increasing amplitude G(t), in order to show the ability of the model in retrieving the dierent nonlinear phenomena experimentally observed. Two congurations are selected for that purpose, U and E0, with a forcing frequency of 115 Hz. This forcing frequency has been arbitrarily chosen to underline how the nonlinear eects appear, it is nonetheless close to 120 Hz which corresponds to an eigenfrequency of conguration E0; it is also close to 109 Hz which corresponds to an eigenfrequency of conguration U. The forcing amplitude G(t) is linearly increasing from 0 to 20 N in 20 s. The spectra of output displacement at the excitation point (x 0 =0.6 m, y 0 =0.01 m) are plotted in Fig. 7. Fig. 7(a) shows that the uniform beam behaves linearly as expected, since the spectrum only contains the excitation frequency. It can be observed in Fig. 7(b) that the spectrum of E0 contains two distinct regimes: there is rst a periodic regime showing an enrichment of the spectrum with odd harmonics, which qualitatively reproduces the results of the experiment of Sec. 2.2. Note that, as compared to the experiments, numerical spectra display only odd harmonics while experimental spectra contain both even and odd harmonics. The presence of only odd harmonics in the simulation is related to the fact that the model contains only cubic nonlinearities. Experimentally, geometric imperfections are unavoidable in real plates, especially in the ABH case where the thickness at the termination is very small. The presence of these imperfections creates in turn quadratic nonlinearity, which gives rise to even harmonics [START_REF] Camier | Non-linear vibrations of imperfect free-edge circular plates and shells[END_REF][START_REF] Ostiguy | Eects of initial geometric imperfections on dynamic behaviour of rectangular plates[END_REF]. The second regime observed in Fig. 7(b) is a chaotic one, associated with the Wave Turbulence (WT) phenomenon, as discussed in Sec. 2.2. Note that the von Kármán model is a moderate perturbation of the motions around the equilibrium position, with linearized angles. Hence the model is not able to produce periodic orbits around other positions than the conguration at rest. Consequently in all the simulated cases, the return to periodic orbits sometimes observed experimentally when still increasing the forcing amplitude, has not been retrieved. The choice of not considering imperfections in the model has been made both for the sake of simplicity and also mainly because in turbulent vibratory regimes with large amplitude, the cubic nonlinearities dominate the quadratic ones, as shown theoretically and numerically in [START_REF] Ducceschi | Dynamics of the wave turbulence spectrum in vibrating plates: A numerical investigation using a conservative nite dierence scheme[END_REF].

(a) (b)

Low frequency ABH performance

In order to observe the potential for an increased eciency of the ABH in low frequency, the beams U, E0 and E0D0 are excited with a ltered white noise on the 100-500 Hz range during 10 s, for several amplitudes, i.e.

several standard deviation σ, as it is done experimentally in Sec. 2.3. The Power Spectral Density of the output velocity of beams E0 and E0D0 at the excitation point (x 0 =0.6 m, y 0 =0.01 m) is displayed in Fig. 8 for dierent amplitudes of excitation. For the covered ABH beam E0D0, experimental results are qualitatively reproduced since energy leaks outside the excitation range, and the transfer becomes more important with the amplitude level, without however aecting the resonance peaks in the excitation range. This energy transfer towards higher frequency (energy cascade due to the Wave Turbulence [START_REF] Touzé | Transition scenario to turbulence in thin vibrating plates[END_REF][START_REF] Auliel | Wave turbulence buildup in a vibrating plate[END_REF]) is much more noticeable and important in the E0 case, since the level outside excitation range easily reaches the level within, and a larger frequency range is concerned (see Fig. 8(c)). At large amplitude, the level of the resonance peaks within the excitation range is slightly reduced because this energy is transferred to higher frequencies where it is dissipated. It is important to note that while very little damping is present in the beam, nonlinearities allow for an interesting reduction of vibration in this case.

In order to evaluate the performance gain of the ABH due to the aforementioned nonlinear eects, an indicator I c,σ is built based on the responses to the white noise excitation of Fig. 8. It is dened as:

I c,σ = 500 100 S c,σ aa (f ) S c,σ pp (f ) df, (37)
where c refers to the current conguration, S aa is the PSD of acceleration output a(t) and S pp is the PSD of the forcing p(t). The comparison is drawn out with a reference case denoted as I ref . This reference case is selected as the uniform beam, conguration U, excited with the minimum amplitude (σ = 0.11N ), hence giving rise to the indicator denoted as E c,σ :

E c,σ = 10 log 10 I c,σ I ref . (38)
Indicator E c,σ is plotted in Fig. 9 for the three congurations U, E0 and E0D0. A quasi-constant value is obtained for the uniform beam, as it is expected. Although the PSD is averaged over 10 s, the variability of the white noise excitation induces that the same exact value cannot be obtained for a linear case at dierent excitation levels, as it is shown by the small variations of case U. The value of the indicator for beam E0 is constant for small σ and then decreases, indicating that less energy is present in the 100-500 Hz range, as it was shown by the spectrum of Fig. 8(c). Finally, the indicator has a constant and smaller value for beam E0D0, reecting a better performance due to the damping layer but also smaller nonlinear eects and no signicant energy transfer towards the high frequency domain, as it is shown in Fig. 8(b). This rst result shows that a compromise has to be found between two competing eects : nonlinearity and damping. A too large value of damping prevents the nonlinear eect to settle down and thus creates an overdamped situation with no dependence on vibration amplitude. The next sections are devoted to a quantitative study of these two eects in order to optimize the low-frequency eciency of the ABH.

Eect of the thickness prole

The importance of the nonlinear eects is rst investigated by simulating the responses of congurations E1 to E4, which have increasing additional length of constant and minimum thickness. Increasing this length should enhance the nonlinear behavior. Moreover, the modal density will also increase in the excited frequency band, hence favoring the turbulent behavior and the energy transfer.

The procedure described in Sec. 4.4 is thus applied to beams E1 to E4 with gradually increasing length of extensions, according to Fig. 4. The results for the indicator E c,σ are plotted in Fig. 10. Conguration E1 shows that adding a small extension positively aects the nonlinear behavior since the structure becomes more sensitive to the excitation level: Indicator E c,σ decreases with respect to the amplitude of excitation and reaches values close to -9 dB. Congurations E2 to E4 are even more interesting from this point of view since the indicator reaches values around -12 dB. However, the similar values obtained for these congurations seem to indicate a saturation eect, i.e. that a longer extension does not result in a better energy transfer or a better dissipation due to nonlinearities.

As known from wave turbulence theory, nite-size eects alters the energy cascade due to the fact that a continuum of wavevectors is not at hand anymore, so that the resonance relationships shall then be veried on a discrete set composed of the modes of the system [START_REF] Cadot | Wave turbulence in vibrating plates[END_REF][START_REF] Kartashova | Discrete wave turbulence[END_REF][START_REF] Miquel | Nonlinear dynamics of exural wave turbulence[END_REF]. In our case where energy transfer has to be optimized, the best situation would then be to have a nonlinear termination as large as possible in order to ensure the best modal density in the low-frequency range, and thus optimize the energy ux toward high frequencies. This situation is out of range for the selected beams where the number of eigenmodes in the tested frequency band 100-500 Hz increases from 8 (E0 conguration) to 17 (conguration E4). Fig. 11 shows the energy spectra for beams E0 to E4, and for the largest excitation amplitude tested, σ=2.6 N, highlighting the fact that the energy is redistributed over the prominent peaks dened by the eigenfrequencies, when extending the additional length, which explains the observed saturation on the indicator. However, better performance could be awaited for two-dimensional ABH as described in [START_REF] Georgiev | Damping of structural vibrations in beams and elliptical plates using the acoustic black hole eect[END_REF][START_REF] O'boy | Damping of exural vibrations in circular plates with tapered central holes[END_REF], where the extension of the minimum thickness region can be signicantly larger.

Balance of damping and nonlinearity

As shown in the previous section, the additional length l add has to be chosen as large as possible in order to enhance the nonlinear eects. The length of the damping layer is now investigated, in order to nd the best compromise between nonlinearity and damping, and to ensure an increasing performance of the ABH thanks to the nonlinear behavior.

In order to observe the competition between damping and nonlinearities in the ABH beam, dierent lengths of damping layer are added to conguration E4, resulting in congurations E4D0 to E4D4 (see Fig. 4). The behavior of the indicator E c,σ for the associated simulations with a ltered white noise excitation are plotted in Fig. 12. Beam E4D0 reveals an interesting combination since the indicator E c,σ is lower than for E4 for small excitation levels, due to the more important damping in the structure. When increasing the vibration amplitude, indicator E c,σ for E4D0 shows then a decreasing behavior indicating that the nonlinearity is activated and a better eciency of the ABH is obtained. The eect can clearly be seen in Fig. 13, which reveals reductions of the PSD down to 14 dB in the 100-500 Hz range (see Fig. 13(b)), with the σ 2 normalization. More specically one can observe that the resonance frequencies tend to larger values when increasing the nonlinearity, reecting the hardening behavior of the plate. Secondly, all resonance peaks widen due to the appearance of the strongly nonlinear behavior associated with the chaotic dynamics, and globally decrease as energy has own to higher frequencies.

When the damping layer length is increased, gradually covering the nonlinear extension, nonlinear eects tend to disappear, as shown in Fig. 12 for E4D2 to E4D4 congurations, where the indicator is seen to stay at a constant value whatever the excitation level. A longer layer, hence a more important damping, however results in a better performance of the ABH. This study shows that both damping and nonlinear eect can be decoupled as the additional length (or the length of the ABH region) and the length (or the material characteristics) of the damping layer can be independently adjusted in order to nd an optimal behavior, depending of the context.

Conclusions and perspectives

This study is devoted to the nonlinear eects in ABH beams, and the potential benet that can be obtained from using them as a vector for transferring energy from the low frequencies, where ABH is known to be usually inecient, to the high frequencies; this transfer is provided by the wave turbulence regime and its associated energy cascade. Experimental results on an ABH beam have been rst reported, evidencing that nonlinearities are activated, giving rise to typical nonlinear phenomena. The thinness of the ABH termination are the support for the geometric nonlinearity, given the amplitude of motions involved.

A von Kármán plate model with variable thickness has then been derived in order to study numerically the main eects involved, with the aim of optimizing the energy transfer and improving the eciency of the ABH in the low frequency range. A full numerical model has been developed, including the geometry variations, the damping due to the added viscoelastic layer, and the geometric nonlinearity due to large-amplitude vibrations.

A modal approach with an energy-conserving scheme has been used for the resolution, this method allowing one to ensure a very good accuracy and convergence of the model for reasonable computational burden, together with a ne description of the frequency-dependent losses present in the structure.

A parametric study on a number of selected geometries has then been led, in order to investigate how the two competing eects represented by both the nonlinearity and the damping, can be tuned for obtaining an important gain in eciency for the ABH. A good compromise has been found by using a long additional termination such that the nonlinearity can develop into a wave turbulence regime, together with a moderate length of the damping layer. Indeed, using a too large value of the damping overdamps all vibrations and counteracts the nonlinear behavior such that no improvement in eciency is obtained with increased amplitudes. On the other hand, letting the chaotic dynamics settle down allows for the energy transfer in high frequencies where vibrations are damped, and thus improves the eciency as amplitude increases. More generally, the results presented in this study show that previous linear models for the ABH may not apply anymore to some particular congurations, or need to be used with great care as they neglect an important nonlinear eect present in the dynamics. Further work dealing with plates embedding two-dimensional circular ABH with nonlinear eects could also be very interesting and closer to possible applications in mechanical engineering. As noted in the course of this study, two-dimensional ABH can present larger extensions with minimal thickness, ensuring a larger modal density in the low-frequency range. This fact should optimize the energy cascade, and could be more advantageous for evacuating the energy toward high frequencies.

Figure 1 :

 1 Figure 1: Spectrograms of output acceleration for the reference beam and the ABH beam excited at various frequencies. (a) reference beam, 83 Hz, (b-d) : ABH beam, excited respectively at 83 Hz (b), 113 Hz (c), and 102 Hz (d). The forcing amplitude is raised from 0 to 15 N in 40 s. Color scale in dB (ref. 1 m.s -2 )

Figure 2 :

 2 Figure 2: Power spectral density of acceleration response to a 100-500 Hz white noise of (a) Uniform beam for σ=0.39 N (yellow), 0.62 N (red), 0.84 N (brown) and 1.4 N (black) and (b) ABH beam for σ=0.26 N (yellow), 0.50 N (red), 0.68 N (brown) and 1.1 N (black), normalized by the variance σ 2 of the excitation, in dB (ref. 1 m 2 .s -4 .Hz -1 .N -2 ). Vertical lines indicate the excitation range. (c) Zoom of (b) on the 100-500 Hz range.

Figure 5 :

 5 Figure 5: (a) transverse modal shape for mode (23,0) of the E0 beam. (c)-(d) in-plane modal shapes for E0, (c) mode Ψ 101,2 belonging to family 1 and (d) mode Ψ 61,3 belonging to family 2.

Figure 6 :

 6 Figure 6: Coupling term Γ p p,p,p for conguration U (black dashed line) and E0 (red full line). (a) Convergence of Γ p p,p,p for p=29 (which correspond to transverse mode (23,0)), as function of the in-plane modes, sorted by eigenfrequencies of increasing order. (b) same convergence, now as function of the in-plane modes sorted by families. Vertical line on (b) separates the rst family of modes Ψ of E0 from the others.

Figure 7 :

 7 Figure 7: Spectrograms of output velocity at point (x 0 =0.6 m, y 0 =0.01 m) for a forcing at 115 Hz from 0 to 20 N in 20 s, for beams (a) U and (b) E0. Color scale in dB (ref. 1 m.s -1 )

Figure 8 :

 8 Figure 8: Power spectral density of velocity response to a 100-500 Hz white noise of (a) U for σ=0.11 N (yellow), 0.98 N (red) and 2.2 N (black), (b) E0D0 for σ=0.11 N (yellow), 1.1 N (red) and 2.4 N (black), and (c) E0 for σ=0.13 N (yellow), 0.99 N (red) and 2.7 N (black), normalized by the variance σ 2 of the excitation, in dB (ref. 1 m 2 .s -2 .Hz -1 .N -2 ).

Figure 9 :

 9 Figure 9: Indicator of performance E c,σ for congurations U (black dashdotted line •), E0 (black dashed line ) , E0D0 (red full line ).

Figure 10 :

 10 Figure 10: Indicator of performance E c,σ for congurations with dierent additional lengths : case U (black dashdotted line •), E0 (black dashed line ), E1 (brown full line ), E2 (red full line ), E3 (yellow full line ) and E4 (black full line ♦).

Figure 11 :

 11 Figure 11: Power spectral density of velocity response to a 100-500 Hz white noise of E0 (blue dashed line), E1 (brown dashdotted line), E2 (red full line), E3 (yellow full line), E4 (black thick full line) for σ=2.6 N, normalized by the variance σ 2 of the excitation, in dB (ref. 1 m 2 .s -2 .Hz -1 .N -2 ).

Figure 12 :

 12 Figure 12: Indicator of performance E c,σ for structures with dierent damping layer length : case U (black dashdotted line •), E4 (black full line♦) E4D0 (red dashdotted line ♦), E4D1 (brown full line ), E4D2 (brown dashed line ), E4D3 (yellow full line ) and E4D4 (yellow dashed line ).

Figure 13 :

 13 Figure 13: (a) Power spectral density of velocity response to a 100-500 Hz white noise for conguration E4D0 for σ=0.12 N (yellow), 0.81 N (red) and 2.6 N (black), normalized by the variance σ 2 of the excitation, in dB (ref. 1 m 2 .s -2 .Hz -1 .N -2 ). (b) Zoom of (a) on the 100-500 Hz range.

Table 2 :

 2 Values of the parameters specic to the congurations simulated with the numerical model.

	U E0 E0D0 E1 E2 E3 E4 E4D0 E4D1 E4D2 E4D3 E4D4 l add [mm] / 0 0 50 100 150 200 200 200 200 200 200 x l1 [mm] / / 0 / / / / 200 150 100 50 0 x l2 [mm] / / 200 / / / / 400 400 400 400 400
	U	E0D0
	E0	E4D0
	E1	E4D1
	E2	E4D2
	E3	E4D3
	E4	E4D4
	Figure 4: Schematic representations of the ABH extremities of the simulated beams. As a memory aid, E stands for extension, D stands for damping and the number immediately after the letter is an indication for the
	concerned length.

A Finite dierence operators

The spatial nite dierence operators δ x-, δ x+ , µ x-, µ x+ and δ yy appearing in Eq. ( 28) are dened as [START_REF] Bilbao | Numerical sound synthesis : nite dierence schemes and simulation in musical acoustics[END_REF]:

In Eq. ( 34), the nite dierence operators in time are dened as:

e t-q n = q n-1 .

(48)


Anaël Grandjean
email: anael.grandjean@lirmm.fr

Victor Poupet
email: victor.poupet@lirmm.fr

L-Convex Polyominoes are Recognizable in Real Time by 2D Cellular Automata

A polyomino is said to be L-convex if any two of its cells are connected by a 4-connected inner path that changes direction at most once. The 2-dimensional language representing such polyominoes has been recently proved to be recognizable by tiling systems by S. Brocchi, A. Frosini, R. Pinzani and S. Rinaldi. In an attempt to compare recognition power of tiling systems and cellular automata, we have proved that this language can be recognized by 2-dimensional cellular automata working on the von Neumann neighborhood in real time. Although the construction uses a characterization of L-convex polyominoes that is similar to the one used for tiling systems, the real time constraint which has no equivalent in terms of tilings requires the use of techniques that are specific to cellular automata.

Introduction

Two-dimensional cellular automata and tiling systems are two different models that can be considered to recognize classes of two-dimensional languages (or picture languages). Although they share some similarities such as locality and uniformity, the two models are fundamentally different.

Tiling systems as language recognizers were introduced by D. Giammarresi and A. Restivo in 1992 [START_REF] Giammarresi | Recognizable picture languages[END_REF] and are based on the model of tile sets introduced by H. Wang [START_REF] Hao | Proving theorems by pattern recognition II[END_REF]. The strength of the model lies in its inherent non-determinism. The system itself is a set of local rules describing valid image patterns and a picture language is recognized by the system if it is the image by a projection of the set of configurations that verify all local rules.

Cellular automata on the contrary are deterministic dynamical models. Introduced in the 1940s by S. Ulam and J. von Neumann [START_REF] John Von Neumann | Theory of Self-Reproducing Automata[END_REF] to study self replication in complex systems they were rapidly considered as computation models and language recognizers [START_REF] Alvy | Real-time language recognition by one-dimensional cellular automata[END_REF]. Contrary to some other classical computation models that inherently work on words, they can be considered naturally in any dimension (the original cellular automata studied by Ulam and von Neumann were 2-dimensional) and are therefore particularly well suited to picture languages. Language recognition is performed by encoding the input in an initial configuration and studying the (deterministic) evolution of the automaton from that configuration. Time and space complexities can be defined in the usual way.

Because tiling systems lack dynamic behavior, some picture languages that can be recognized by cellular automata with minimal space and time complexity (in real time) cannot be recognized by tiling systems, such as the language of square pictures with vertical symmetry.

Conversely, the non-determinism of tiling systems should allow the recognition of languages that cannot be recognized by cellular automata in low time complexities. It is straightforward for instance to verify that the language considered in [START_REF] Terrier | Two-dimensional cellular automata recognizer[END_REF] as an example of language that cannot be recognized in real time by a cellular automaton working on the Moore neighborhood but can be recognized on the von Neumann neighborhood can be recognized by a tiling system, thus proving that tiling systems and real time cellular automata on the Moore neighborhood are incomparable.

Because the language of L-convex polyominoes was recently proved to be recognizable by tiling systems when it was previously though not to be, we decided to investigate its recognizability by real time von Neumann neighborhood cellular automata. Although the language was also recognized by cellular automata, the construction turned out to be quite different from the case of tiling systems and used some techniques specific to cellular automata (and possibly von Neumann neighborhood cellular automata). This article describes said construction.

Definition 2 (Configuration

). A d-dimensional configuration C over the set of states Q is a mapping from Z d to Q.
The elements of Z d will be referred to as cells and the set of all d-dimensional configurations over Q will be denoted as Conf d (Q).

Given a CA A = (d, Q, N , δ), a configuration C ∈ Conf d (Q) and a cell c ∈ Z d , we denote by N C (c) the neighborhood of c in C : N C (c) : N → Q n → C(c + n)
∆(C) = Z d → Q c → δ(N C (c))
The action of the global transition rule makes A a dynamical system over the set Conf d (Q). Because of this dynamic, in the following we will identify the CA A with its global rule so that A(C) is the image of a configuration C by the action of the CA A, and more generally A t (C) is the configuration resulting from applying t times the global rule of the automaton from the initial configuration C.

Picture Recognition

From now on we will only consider 2-dimensional cellular automata (2DCA), and the set of cells will always be Z 2 .

Definition 4 (Picture). For n, m ∈ N and Σ a finite alphabet, an (n, m)picture (picture of width n and height m) over Σ is a mapping

p : 0, n -1 × 0, m -1 → Σ Σ n,m
denotes the set of all (n, m)-pictures over Σ and Σ * , * = n,m∈N Σ n,m the set of all pictures over Σ. A picture language over Σ is a set of pictures over Σ.

Definition 5 (Picture Configuration). Given an (n, m)-picture p over Σ, we define the picture configuration associated to p with quiescent state q 0 / ∈ Σ as

C p,q0 :    Z 2 → Σ ∪ {q 0 } x, y → p(x, y) if (x, y) ∈ 0, n -1 × 0, m -1 q 0 otherwise
Definition 6 (Picture Recognizer). Given a picture language L over an alphabet Σ, we say that a 2DCA A = (2, Q, N , δ) such that Σ ⊆ Q recognizes L with quiescent state q 0 ∈ Q \ Σ and accepting states Q a ⊆ Q in time τ : N 2 → N if, for any picture p (of size n × m), starting from the picture configuration C p,q0 at time 0, the origin cell of the automaton is in an accepting state at time τ (n, m) if and only if p ∈ L. Formally,

∀n, m ∈ N, ∀p ∈ Σ n,m , A τ (n,m) (C p,q0 )(0, 0) ∈ Q a ⇔ p ∈ L
Because cellular automata work with a finite neighborhood, the state of the origin cell at time t (after t actions of the global rule) only depends on the initial states on the cells in N t , where N 0 = {0} and for all n, N n+1 = {x + y, x ∈ N n , y ∈ N }. The real time function is informally defined as the smallest time such that the state of the origin may depend on all letters of the input :

Definition 7 (Real Time). Given a neighborhood N ⊂ Z d in d dimensions, the real time function τ N : N d → N associated to N is defined as τ N (n 1 , n 2 , . . . , n d ) = min{t, 0, n 1 -1 × 0, n 2 -1 × . . . × 0, n d -1 ⊆ N t }
When considering the specific case of the 2-dimensional von Neumann neighborhood, the real time is defined by τ NvN (n, m) = n + m -2. There is however a well known constant speed-up result : Proposition 1 (folklore). For any k ∈ N, any language that can be recognized in time (τ NvN +k) by a 2DCA working on the von Neumann neighborhood can also be recognized in real time by a 2DCA working on the von Neumann neighborhood. So it will be enough to prove that a language is recognized in time (n, m) → n + m + k for some constant k to prove that it is recognized in real time.

Polyominoes

Definition 8 (Polyomino). A placed polyomino is a finite and 4-connected subset of Z 2 . A polyomino is the equivalence class of a placed polyomino up to translation.

Definition 9 (HV-Convexity). A polyomino p is said to be horizontally (resp. vertically) convex if any cell between two cells of the polyomino on a same horizontal (resp. vertical) line is also a cell of the polyomino :

∀x 1 , x 2 , x 3 , y ∈ Z, x 1 ≤ x 2 ≤ x 3 ∧ (x 1 , y) ∈ p ∧ (x 3 , y) ∈ p ⇒ (x 2 , y) ∈ p
A polyomino is HV-convex if it is both horizontally and vertically convex (see Figure 2).

We will now present the notion of L-convex polyomino, first introduced in [START_REF] Castiglione | Reconstruction of l-convex polyominoes[END_REF] to classify HV-convex polyominoes. Informally, an L-convex polyomino p is such that for any two of its cells there exists a 4-connected path of cells of p that connects them such that the path changes direction at most once (see Figure 3).

The following remarks will lead to a formal definition of L-convex polyominoes : Fig. 3. The polyomino on the left is Lconvex (the figure shows an inner path connecting two cells with at most one direction change, and there is such a path for any pair of cells). The polyomino on the right is HVconvex but not L-convex as illustrated by the pair of highlighted cells for which there is no inner connecting path that changes direction at most once.

0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Fig. 4. A polyomino (left) and its corresponding picture over {0, 1} (right). When this picture is encoded as a configuration of a cellular automaton, the origin of the automaton is on the lower left corner of the picture.

a path that changes direction at most once connecting two cells of a polyomino p on the same row (resp. column) is fully horizontal (resp. vertical) therefore L-convex polyominoes are HV-convex ;

-if c 1 = (x 1 , y 1 ) and c 2 = (x 2 , y 2 ) are two cells in an L-convex polyomino p, either a 1 = (x 1 , y 2 ) or a 2 = (x 2 , y 1
) is a cell of p because a 1 and a 2 are the angles of the only two paths connecting c 1 and c 2 that change direction at most once ; if a polyomino p is HV-convex and such that for any two of its cells c 1 = (x 1 , y 1 ) and c 2 = (x 2 , y 2 ) either a 1 = (x 1 , y 2 ) or a 2 = (x 2 , y 1 ) is a cell of p, then p is L-convex since by HV-convexity, the whole path connecting c 1 to c 2 going through a 1 or a 2 is in p.

Definition 10 (L-Convexity). A polyomino p is L-convex if it is HV-convex and verifies

∀x 1 , x 2 , y 1 , y 2 ∈ Z, (x 1 , y 1 ) ∈ p ∧ (x 2 , y 2 ) ∈ p ⇒ (x 1 , y 2 ) ∈ p ∨ (x 2 , y 1 ) ∈ p
Given a polyomino p, the picture over the alphabet {0, 1} associated to p is the picture whose dimensions are the dimensions of the minimal bounding rectangle of p, where the cell has state 1 if the corresponding cell is in the polyomino and 0 otherwise (see Figure 4). We define the language L L-convex as the picture language of all L-convex polyomino pictures.

Main Result

This section will be entirely devoted to the proof of the following result Theorem 1. The picture language L L-convex of L-convex polyomino pictures is recognizable in real time by a 2DCA working on the von Neumann neighborhood.

The proof will be done by describing the behavior of a 2DCA working on the von Neumann neighborhood that recognizes L L-convex in real time. In this description we will use cardinal directions north, south, east and west to denote the different directions on the configuration as follows :

north is towards the increasing y axis ; south is towards the decreasing y axis ; east is towards the increasing x axis ; west is towards the decreasing x axis.

With such conventions, the origin of the automaton is located at the south-west (SW) angle of the picture in the initial configuration and the picture therefore extends from the origin eastward and northward.

Preliminary Check

First of all, the automaton must check that the input is the picture of a HVconvex polyomino.

To do so, during the first step of the computation, each cell containing a 1 considers its neighbors and remembers which of them also contains a 1. Then a signal moves westward from the eastmost point of each row and southward from the northmost point on each column. These signals check that each row and each column contains exactly one segment of connected 1 symbols. Moreover, the signals check that the segment of 1 on each line and column is connected to that of the neighbor rows and columns using the neighboring information gathered during the first step.

These two properties guarantee that the polyomino is connected, HV-convex and that the picture's dimensions are that of the minimal bounding rectangle (no empty row or column). If an error is found on a row or column, the signal is directed towards the origin and the input is not accepted.

We can now assume that the input corresponds to a HV-convex polyomino picture, and must determine whether it is also L-convex.

Characterization of L-Convex Polyominoes

We will now present the characterization of L-convex polyominoes that will be used by the automaton. It is a slighly rephrased version of the characterization presented in [START_REF] Brocchi | A tiling system for the class of L-convex polyominoes[END_REF] (Theorem 2).

Given a polyomino p, we say that a cell of p is a corner if it has two consecutive neighbors that are not in p. We classify corners depending on the directions in which such neighbors not in p are located : a north-east (NE) corner is one such that the northern and eastern neighbors are not in p, and we similarly have NW , SW and SE corners (see Figure 5 for an illustration of NE corners). Note that corner types are not exclusive : a cell can for instance be both a NE and NW corner.

Proposition 2 (Characterization of L-convex polyominoes [START_REF] Brocchi | A tiling system for the class of L-convex polyominoes[END_REF]). A HVconvex polyomino p is L-convex if and only if for every NE corner c = (x, y), denote by (x, y ) the southest cell of p on the same column as c, and (x , y) the westmost cell of p on the same row as c, there is no cell (x , y ) of p verifying any of the following three conditions (a). x > x (resp. x > x ) and y < y (b). x < x (resp. x < x ) and y > y (c). x < x (resp. x < x ) and y < y and the symmetric conditions holds for all NW corners (in the South and East directions).

NE corners

(x, y) (x', y) (x, y') (x', y') a b c h 1 v 1 h 2 v 2
Fig. 5. A HV-convex polyomino is L-convex is for any of its NE corners (represented as dark grey cells), no cell of the polyomino lies in any of the three zones represented in hatched light grey, and symmetrically for all of its NW corners. The illustrated polyomino is not L-convex because there are two cells in the lower left hatched area (these cells cannot be connected to the represented NE corner by an inner path with at most one direction change).

Proof (sketch). It is enough to verify that all pairs of corners of a HV-convex polyomino are connected by an inner path with at most one change of direction. Moreover by symmetry we can consider only NW and NE corners.

The cells (x , y), (x, y ) and (x , y ) in the characterization represent the farthest points that can be reached from a given corner in their respective directions. Cells of the three restricted areas cannot be connected to the corner and conversely all cells not in these areas can be connected to the corner.

Note that because the polyomino is assumed to be HV-convex it is enough to check that there is no polyomino cell on the two lines extending from the starting check point (represented in dark hatched grey in Figure 5). For instance, for the condition (a), it is enough to check that there is no cell (x , y -1) with x > x and no cell (x + 1, y ) with y < y in the polyomino. This follows from the 4-connectedness of the polyomino.

Although the conditions to verify are perfectly symmetric for NE and NW corners, when implementing it on a real time cellular automaton the case of NE corners is significantly simpler because all signals move towards the origin at maximum speed so the result of the verification easily arrives on time. On the other hand, for NW corners, some signals move eastward (away from the origin) so it would take too much time to send the signal all the way to the east side and back to the origin. We will therefore now focus on implementing the characterization for NE corners and come back to the NW corners at the end of the proof.

Compression and Marking

The characterization from Proposition 2 depends on cells being able to tell if there is a polyomino cell in a given direction from them. To make sure that each cell knows this information, consider signals going eastward from the west side of each row. If the initial configuration is the picture configuration of a HV-convex polyomino, there is exactly one segment of 1 symbols on each row. Before the signal meets the first 1, cells can be notified that there is no 1 westward and that there is at least one 1 eastward. On the segment of 1, cells are notified of whether they are a border cell or an inner cell and, after the segment of 1, all cells are notified that there is a 1 westward and none eastward. Of course the same thing can be done on columns with northward signals.

Now consider a horizontal compression of the input as illustrated by Figure 6. To compress the input, consider that each cell can now hold two initial states instead of one (this can be done by increasing the number of states of the automaton) and move all states westward unless the column in which they should go is full (contains two states) or is out of the boudaries of the initial configuration (ignore the darker dots from Figure 6 for the moment).

Such a compression takes n 2 time steps where n is the width of the input. During these steps, no signal can propagate westward because the initial data is already moving west at maximum speed but the time lost performing the compression can be recovered afterwards because each cell now sees twice as many states horizontally, which means that relatively to the original states, horizontal signals can perform two steps at a time.

During the compression, signals can however be propagated eastward (as illustrated by the darker dots in Figure 6). This means that while the compression is taking place, the signal indicating to each cell if it has 1 symbols east or west can propagate, so that at the end of the compression, cells have access to this information.

By performing a vertical compression after the horizontal one we can otain in half of the real time a compressed copy of the initial configuration on which every cell now has the added information of whether there is a 1 in any of the t=0 t=1 t=4

... four directions. Moreover, later in the construction we will need to know which columns correspond to the same horizontal segment on the southern border of the polyomino, so we also propagate northward signals during the vertical compression from the borders of all horizontal segments of the southern border of the polyomino (see dashed northward arrows in Figure 8). After both compressions, the computation of the automaton can properly start and in this computation horizontal and vertical signals can propagate twice as fast and all information is twice closer to the origin. This means that the compressed run of the automaton can behave exactly as if the configuration was not compressed but was given the extra information propagated by the eastward and northward signals from the beginning 1 . We will now ignore the compression in the following explanations, and simply consider that the information propagated by the northward and eastward signals is readily available to each cell.

Remark: As it is described, it looks as if cells should know when the compression is finished to start performing the next task (be it the second compression or the accelerated simulation of the uncompressed automaton). However, one can show that cells can asynchronously start the next task as soon as they have the necessary information to do so. It is sufficient to detect when all cells in their neighborhood have finished the compression to perform one step of the next task. From there, we can show that if each cell advances the following task as soon as it has enough information to do so, cells that have completed the compression early will be slowed down progressively to wait for the further cells to catch up. However, the last cells to finish the compression will never be slowed down as all other cells have the necessary information available to them. This means that by continuing the computation after the compression as soon as the information is available, all cells are at least as advanced as if all had started their computation at the time when the compression is finished, thus negating the need to synchronize all cells after the compression. 1 This compression technique works in our case because the automaton (as it will be described later) only uses horizontal and vertical signals that change directions a bounded number of times. It is only possible to simulate two steps of the uncompressed automaton if they only involve horizontal or vertical movement, not both.

First Conditions of the Characterization

With the informations we have, checking conditions (a) and (b) of Proposition 2 is very easy as it is only a matter of sending a westward and a southward signal from each NE corner. When these signals reach the border of the polyomino, they check that there are no 1 in the corresponding area by using the information that was transmitted to each cell during the compressions. If a 1 is found where it should not be, a signal is directed towards the origin to indicate that the input should not be accepted.

There are no conflicting signals during this step because there can be at most one NE corner per column and one at most per row.

The Third Condition

The third condition from Proposition 2 is much more complex to implement. It requires sending a westward signal h 1 and a southward signal v 1 from each NE corner and having these signals generate secondary signals h 2 (westward, from the collision of v 1 with the border of the polyomino) and and v 2 (southward from the collision of h 1 and the border). The intersection of h 2 and v 2 indicate the cell on which condition (c) should be checked, as illustrated by Figure 5.

Two problems arise when implementing this behavior :

although v 1 and h 1 signals originating from different NE corners will never overlap, if two signals arrive on the same row or column they will produce v 2 or h 2 signals that might overlap ; h 2 signals might intersect with many v 2 signals, but only one of them originates from the same NE corner. It is therefore necessary to ensure that the verification of condition (c) is not performed on cells that do not correspond to a valid intersection of h 2 and v 2 signals.

Priority Rule To solve the first problem, we use a simple priority rule : if two NE corners c 1 and c 2 are north of the same horizontal segment on the southern border of the polyomino, we can ignore the easternmost one. There are two cases to consider (illustrated by Figure 7). Assume c 1 lies north-west of c 2 :

if the westward h 1 signal from c 1 reaches the border east of that from c 2 (left of Figure 7), then the area that would be checked by considering the intersection of the signals v 2 and h 2 from c 1 (dark grey area in the Figure )  is east of the one that would be considered by the intersection from c 2 (light grey area) and therefore contains it entirely, which means that it is not necessary to check the area indicated from c 2 ; if on the contrary the h 1 signal from c 1 arrives west of that from c 2 (right of Figure 7), the HV-convexity of the polyomino guarantees that there can be no 1 in either of the two areas considered by the intersections from c 1 and c 2 since there is at least one 1 north west of where the h 1 signal from c 2 arrives, no 1 west of that point so there cannot be any 1 west and south of it. In this case, it doesn't matter which intersection is considered since neither will find a contradiction with the (c) condition from Proposition 2. A symmetrical argument shows that it is sufficient to consider signals originating from the southernmost of two NE corners whose h 1 signals arrive on the same vertical segment of the western border of the polyomino. Horizontal and vertical signals are however handled differently because the last part of the construction is not symmetrical.

We want to make sure that there are as many v 1 signals as there are distinct (non-overlapping) h 2 signals. To do so, v 1 signals are not sent directly by NE corners but rather sent by the h 1 signal when the h 1 signal knows that the corner it originated from is the westmost of the corresponding horizontal segment in the southern border (see Figure 8). When an h 1 signal finds a cell of the polyomino north before reaching the border of the southern segment (dashed line in the figure), it knows there is another NE corner west for that segment and therefore disappears. On the contrary, if such a signal reaches the border of the southern segment it sends the v 1 signal southward.

Counters For v 2 signals, we need to solve the second problem that was described previously which is to determine which of the possibly many h 2 signals intersected is the one that originated from the same NE corner. To do so, h 1 signals produced by NE corners will count how many v 2 signals they cross while going west. If a h 1 signal crossed n v 2 signals, then the v 2 signal it produces will consider that its corresponding h 2 signal is the (n + 1)-th to last one (the last n are not the one that should be considered).

Figure 9 illustrates why the result of such a behavior is correct. Consider an NE corner c 2 such that its h 1 signal crossed the is the first that the v 2 signal from c 2 encounters, and the later one should be ignored (lower circled intersection) ; if on the contrary the v 1 signal from c 2 arrives south of that of c 1 (right part of the figure), the real intersection that should be considered is the last one but by considering the first the automaton will not find any contradiction to condition (c) since by HV-convexity of the polyomino there are no 1 south and west of either of the two intersections (so it will pick the wrong intersection but that will not change the final result). In order to implement this rule, signals need to carry a binary counter. This counter should follow the signal at maximal speed, and will be incremented by the h 1 signal for each v 1 signal encountered (which can be easily done as incremental binary counters can be implemented on one-way one dimensional CA). For technical reasons, the counter has an initial value of 1.

As for the v 2 signal, as it crosses h 2 signals it checks if the area south-west of said intersection contains a 1 (which can be done instantly because of the information gathered during the initial compressions), and if so decrements the counter 2 . If the counter is equal to 0 (decreasing a 0 counter leaves it at 0) when the v 2 signal reaches the southernmost border of the picture, no error is detected, but if it is positive then a message is sent to the origin to indicate that the input is not L-convex.

This works because we know that if there is no 1 south-west of a cell c, there is no 1 south-west of a cell south of c. If the h 1 signal crosses n v 1 signals the counter indicates (n + 1) when the v 2 signal starts moving south and if the counter is at 0 it means that the (n + 1) last intersections were correct according to condition (c) and therefore the (n + 1)-th to last was correct.

Checking that the counter is 0 takes log(n) steps where n is the maximal value of the counter (log(n) is the maximal length of the counter). If the counter is incremented to n it means that there are at least n NE corners north of the one from which the signal originated. This means that if the signal moves from this corner towards the origin (south or west) at maximal speed, it would reach the origin at least n steps before the real time, and therefore it can spend log(n) steps checking the value of the counter and still arrive in real time.

Moreover, conflicts of overlapping counters can be resolved by the priority rule described previously. Precedence must always be given to the counter corresponding to the southernmost NE corner :

when an h 1 signal reaches the west border of the polyomino, it marks the cell on which the v 2 signal is produced ; if a v 2 signal moves through such a marked cell, it is erased ; if the counter following a v 2 signal is on a cell where a new v 2 signal is created, the counter is invalidated (the end symbol is erased) so that the new v 2 signal has precedence over it. An invalidated counter will ignore decrementations and will ignore the test to 0 at the end.

The North-West Corners

The previous subsections describe how NE corners can properly implement the characterization of L-convex polyominoes from Proposition 2. NW corners will behave in a very similar fashion, but special care must be taken to prove that the result of their verification can reach the origin in real time. On a regular configuration, a signal issued from a NW corner needs to go east through most of the polyomino, then south and then the result of the verification should travel back west to the origin. In doing so the signal goes twice through the width of the input which cannot be done in real time. To solve 2 Decrementation can also be performed on a binary counter moving at maximal speed but in that case the length of the counter is not reduced when going from an 2 n to (2 n -1), but instead a leading 0 is added, which is not a problem for our construction.

this problem, we consider the path of the signal during a horizontal compression of the configuration :

the signal starts from the NW corner on the cell c = (x, y) ; during the x 2 first steps the signal moves west with the compression, and when the cell is compressed, the h 1 signal is sent eastward ; meanwhile, the cell (x , y) that should have been the target of the h 1 signal moves left with the compression. The h 1 signal and the cell arrive at the cell ( x 2 , y) at time ( x 2 , y) ; the signal v 2 from c moves south until it reaches the southern border of the input after y steps. At this point a delay of at most log(h -y) steps is incurred to check the value of the counter (h is the total height of the input) ; the result of the verification is directed towards the origin, it arrives at time x 2 + y + log(h -y) + x 2 < x + h which is before real time. NW corners can therefore properly perform the necessary verifications to implement the characterization from Proposition 2, which concludes the proof of Theorem 1.

1 (

 1 Cellular Automaton). A cellular automaton (CA) is a quadruple A = (d, Q, N , δ) where d ∈ N is the dimension of the automaton ; -Q is a finite set whose elements are called states ; -N is a finite subset of Z d called neighborhood of the automaton ; δ : Q N → Q is the local transition function of the automaton.

From the local transition

  function δ of a CA A = (d, Q, N , δ), we can define the global transition function of the automaton ∆ : Conf d (Q) → Conf d (Q) obtained by applying the local rule on all cells :

Definition 3 (

 3 Von Neumann and Moore Neighborhoods). In d dimensions, the most commonly considered neighborhoods are the von Neumann neighborhood N vN = {c ∈ Z d , ||c|| 1 ≤ 1} and the Moore neighborhoodN M = {c ∈ Z d , ||c|| ∞ ≤ 1}.Figure 1 illustrates these two neighborhoods in 2 dimensions.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The von Neumann (left) and Moore (right) neighborhoods in 2 dimensions. Fig. 2. Three polyominoes. The center and right ones are vertically convex, the right one is HVconvex.

Figure 5

 5 Figure 5 illustrates this characterization.

Fig. 6 .

 6 Fig.[START_REF] John Von Neumann | Theory of Self-Reproducing Automata[END_REF]. Horizontal compression of the input, with eastward transmission of information (dark dots).

Fig. 7 .

 7 Fig. 7. When two v1 signals arrive on the same row we can always safely ignore the one originating from the eastmost NE corner.

2 Fig. 8 .

 28 Fig. 8. h1 signals from an NE corner are interrupted if they detect that there is another NE corner west whose v1 signal would arrive on the same horizontal segment on the southern border of the polyomino. v1 signals are sent by h1 signals on the westmost column corresponding to the southern horizontal segment.

1 Fig. 9 .

 19 Fig.9. Considering that for each v2 signal crossed by the h1 signal from c2, one of the last intersections with an h2 signal should be ignored by the subsequent v2 signal leads to a correct characterization.