Sébastien Héron

Fabrice Pardo

Patrick Bouchon
email: patrick.bouchon@onera.fr

Jean-Luc Pelouard

Riad Haïdar

Modal method for second harmonic generation in nanostructures

A modal method for the study of second harmonic (SH) generation in nanostructures is depicted both for plane waves and focused beam illumination. It is based on the B-spline modal method, and relies on the undepleted pump approximation. The nonlinear polarization induced by an incident plane wave or a focused beam generates a source term at the doubled frequency. The latter is divided into a finite number of sub-sources, and the SH field is subsequently computed by integration of these sub-source contributions.

INTRODUCTION

In the near-and the mid-infrared, nanostructured devices are giving promising results for controlling the localization of the electric field and for strong intensity field enhancement [START_REF] Bouchon | Total funneling of light in high aspect ratio plasmonic nanoresonators[END_REF][START_REF] Grosjean | Diabolo nanoantenna for enhancing and confining the magnetic optical field[END_REF][START_REF] Novotny | Antennas for light[END_REF][START_REF] Chevalier | Giant field enhancement in electromagnetic Helmholtz nanoantenna[END_REF]. Since nonlinear optics processes benefit from a large enhancement of the local electric field, such nanostructures are very promising [START_REF] Kim | Highharmonic generation by resonant plasmon field enhancement[END_REF][START_REF] Schuller | Plasmonics for extreme light concentration and manipulation[END_REF][START_REF] Genevet | Large enhancement of nonlinear optical phenomena by plasmonic nanocavity gratings[END_REF][START_REF] Kauranen | Nonlinear plasmonics[END_REF]. Their design and optimization need fast numerical simulation tools. So far, various numerical methods have been presented to account for nonlinear processes in periodic nanostructures [START_REF] Bravo-Abad | Modeling nonlinear optical phenomena in nanophotonics[END_REF]. In particular, methods have been developed for computing the second harmonic (SH) generation based on the Green function method [START_REF] Borisov | Second harmonic generation from arrays of sub-wavelength cylinders[END_REF], the Dirichlet-to-Neumann map [START_REF] Yuan | Analyzing second harmonic generation from arrays of cylinders using dirichlet-to-neumann maps[END_REF], the boundary element method [START_REF] Mäkitalo | Boundary element method for surface nonlinear optics of nanoparticles[END_REF], the surface integral method, [START_REF] Butet | Secondharmonic generation from periodic arrays of arbitrary shape plasmonic nanostructures: a surface integral approach[END_REF] or the finite difference time-domain method [START_REF] Francés | Tensorial split-field finitedifference time-domain approach for second-and thirdorder nonlinear materials[END_REF]. In addition, modal methods are known to be efficient for the study of a stack of nanostructured layers, and the Fourier modal method has been adapted for the study of SH generation [START_REF] Nakagawa | Analysis of enhanced second-harmonic generation in periodic nanostructures using modified rigorous coupled-wave analysis in the undepleted-pump approximation[END_REF][START_REF] Bai | Fourier modal method for the analysis of second-harmonic generation in two-dimensionally periodic structures containing anisotropic materials[END_REF][START_REF] Paul | A numerical approach for analyzing higher harmonic generation in multilayer nanostructures[END_REF].

More recently, the B-spline modal method (BMM) was introduced to solve linear electromagnetic problems [START_REF] Edee | Analysis of defect in extreme uv lithography mask using a modal method based on nodal b-spline expansion[END_REF][START_REF] Bouchon | Fast modal method for subwavelength gratings based on b-spline formulation[END_REF][START_REF] Walz | B-spline modal method: a polynomial approach compared to the Fourier modal method[END_REF][START_REF] Granet | Efficient implementation of b-spline modal method for lamellar gratings[END_REF]. In this method, the eigenvalue equation describing the electromagnetic problem is sparse, allowing the calculation speed to increase by solving the problem for only a few eigenmodes. Additionally, the B-splines can be nonuniform which allows for a better description of permittivity breaks. In this article, we adapt the BMM to map the SH field throughout periodic nanostructured layers for plane waves (of arbitrary incident angle) under the undepleted pump approximation. In addition, to simulate experimental setups with a focused incident light, we develop a computation of SH field for a cylindrical focused beam at the fundamental frequency. In Section 2, the consecutive steps involved in the computation are depicted, both for plane waves and focused regimes. Then, in Section 3, numerical examples are described and the convergence of the method is investigated.

THEORETICAL APPROACH

In this section, a layered nanostructure, invariant along the y direction and periodic along the x direction (period d), is considered (see Fig. 1). The layer B, structured along the x direction, is surrounded by two semi-infinite media, A and C. The interface between A and B (resp. B and C) is noted α (resp. β). The incident plane wave is TE-polarized (i.e., the electric field is along the y direction) and impinging on the grating at oblique incidence, making an angle θ with the z axis. The SH wave is also TE-polarized. The theoretical method to compute the generated SH field in periodic nanostructures is detailed first for a single plane wave of frequency ω, and is then generalized to a focused beam.

A. Method for Plane Wave SH Calculation 1. Construction of the Fundamental Field

In layer B, the field at the fundamental frequency is computed thanks to the BMM [START_REF] Bouchon | Fast modal method for subwavelength gratings based on b-spline formulation[END_REF]. In this method, the electric field is interpolated on N samples on the x axis with a nonuniform basis of B-splines, which are piecewise continuous polynomials of degree p [START_REF] Boor | A Practical Guide to Splines[END_REF]. This field is determined through a modal equation M B e m k 2 m e m in TE-polarization, where M B is a N × N matrix deduced from the Maxwell and constitutive equations, the scalar k m is the mth mode wave vector along z, and the vector e m is the complex amplitude of the electric field of the mth mode at the sampling points. The eigenvectors are gathered in a matrix E ω B , whose mth column corresponds to the eigenvector e m . E ω B being computed, H ω B for the x component of the magnetic field can also be calculated. The complex y component E of the electric field and x component ℋ of the magnetic field, at the sampling coordinates x and fixed height z, can be represented as superpositions of two waves, one upgoing and the other downgoing in the structure:

EzE zE -z;
(1)

ℋzℋ z -ℋ -z: (2)
These terms are vectors of size N because of the sampling along x over one period. The fields are written at the time frequency ω, thanks to their projections on the eigenmodes basis:

E zE ω B diage ik m z-z 0 U z 0 ; (3)
ℋ zH ω B diage ik m z-z 0 U z 0 : (4)
To find the total fields, it suffices to find the two vectors of the mode amplitudes U , which are given relatively to a height z 0 , and to propagate them in the eigenmodes basis (e -iωt is omitted). To this end, the scattering matrices at each interface are determined thanks to the continuity of the tangential electric and magnetic fields (see Fig. 2) [START_REF] Li | Multilayer modal method for diffraction gratings of arbitrary profile, depth, and permittivity[END_REF]. For instance, S α BA stands for the transmission scattering matrix from medium A to layer B, and S α BB for the reflection scattering matrix on interface α in layer B.

The electric field is thus obtained in layer B at any z:

E B zE ω B P α z P β z S β BB P α β ΓS α BA U -z α ; (5)
where Γ I -S α BB P β α S β BB P α β -1 , and P α z (resp. P β z ) is the propagator from interface α (resp. β) to height z:P z 1 z 2 diage ik z jz 2 -z 1 j . U -z α is the vector for the incident downgoing wave at interface α in medium A.

Construction of the SH Field

In the following, all the equations and fields are given at the doubled frequency 2ω. The grating layer B comprises a nonlinear material for which it is assumed a second order nonlinear susceptibility of the type χ 2 yyy χ 2 , and χ 2 xyy χ 2 zyy 0. Thus, it induces the complex second order nonlinear polarization [START_REF] Boyd | Nonlinear Optics[END_REF]:

P 2 2ω;zε 0 χ 2 2ω;zE 2 B ω;z: (6)
This term generates a SH field E 2ω , which is calculated in the following. The total field at 2ω is obtained by linear superposition of the fields generated in the structure by independent polarization layers. The undepleted pump approximation is used to state that the fundamental field stays unperturbed through the nonlinear process.

At every height z, a sub-layer of thickness dz acts as a source creating a fraction of the SH field. Thus, to determine the field generated by a single sub-layer at height z, we first consider that it is the only source in the whole system. For the generated fraction of the x component of the magnetic field dℋ at 2ω, the integration of the Maxwell-Ampere equation over regions of thickness dz 0 on a period d gives dℋz -dℋz --2iωP 2 zdz;

where z z dz 2 . Then, Eq. ( 7) is projected on the eigenmode basis thanks to Eqs. ( 2) and ( 4):

dU z -dU -z -dU z -dU -z - -2iωℋ 2ω B -1 P 2 zdz: (8)
In addition, the continuity of the electric field across the sub-layer, and the reflection scatterings at interfaces β and α give the following equations (see Fig. 3): dU z dU -z dU z -dU -z -;

(9) dU z -M β dU -z -; (10)
dU -z M α dU z ; [START_REF] Yuan | Analyzing second harmonic generation from arrays of cylinders using dirichlet-to-neumann maps[END_REF] where M i P i z S i BB P z i are functions of z i α; β. Equations ( 8)-( 11) can be reduced to the following system:

M β dU -z -dΠzdU z ; (12)
M α dU z dΠzdU -z -; (13)
where dΠz-iωP 2 zdz. The contribution of a single sublayer is then given by

dU z I -M β M α -1 M β IdΠz; (14) dU -z -I -M α M β -1 M α IdΠz: (15)
Working in the limit dz → 0, one gets z z -z. Each sub-layer is considered linearly independent from the others. Then, the contributions of all the sub-layers are combined to obtain the value of the total SH field at 2ω:

E zE 2ω B Z z z β P z 0 z dU z 0 ; (16)
E -zE 2ω B Z z α z P z 0 z dU -z 0 : (17)
Nevertheless, the problem must be discretized for numerical processing, preventing the different sub-layers from being infinitely small. The infinitesimal quantity dz must then be replaced by a finite one, Δz.

B. Method for SH in Focused Regime

To design nanostructures, it is more significant to use focused beams as inputs rather than plane waves to take their experimental characterization into account. Here, a method is depicted to numerically model beams focused along the x direction as they can be decomposed in a plane wave basis.

Construction of the Fundamental Field

The structure is supposed d-periodic along x, and the wavelength of the beam is fixed at λ. The incident beam is described by its electric field distribution E i x; z in the xOz plane as it is considered infinite along the y direction (see Fig. 4). To gain access to its plane wave spectrum, it can be written as an inverse Fourier in the reference plane z z 0 :

E i x; z z 0 1 2π Z R Ẽi k x e ik x x dk x ; (18)
where the wave vector is k k x ;k z such that jkj2π∕λ.

The Fourier coefficients are given by

Ẽi k x Z R E i x; z z 0 e -ik x x dx: (19)
An incident beam optically focused on a nanostructure does necessarily contain only propagative orders in its Fourier decomposition. It then reduces the space of integration to k x respecting the condition jk x j ≤ 2π∕λ, ensuring the existence of an associated real k z :

E i x; z z 0 Z jk x j≤2π∕λ Ẽi k x e ik x x dk x : (20)
For the sake of numerical simulations, a discrete spectrum of k x is needed. To this end, a super-period of size D Md is considered. It contains a finite number M of periods, such that D is much larger than the spatial extent of the beam at z 0. This super-period D stands for the repetition period of the focused field: Ex D; zEx; z. This way, the fields can be written as a Fourier series over a discrete spectrum. Equation [START_REF] Edee | Analysis of defect in extreme uv lithography mask using a modal method based on nodal b-spline expansion[END_REF] for the incident D-periodic field E i D takes the form

E i D x; z z 0 X jmj<D∕λ Ẽi m e ik m x x ; (21)
where k m x 2πm∕D and the wave vector is k k m x ;k z such that jkj2π∕λ. Note that the condition for propagating modes now turns to jmj <D∕λ. The Fourier coefficients are given by

Ẽi m 1 D Z D∕2 -D∕2 E i D x; z z 0 e -ik m x x dx: (22)
The fundamental m components of the field E f k m

x x; z, resulting from the interaction between the incoming plane waves e ik m x x , and the nanostructures are then constructed through the modal method previously described. Each mode of the field respects the condition of pseudo-periodicity:

E f m x d; zE f m x; ze ik m x d
. This gives the relative phase between the various mode wave vectors k m

x from one period of the super-period to another. In a given period n, the total field can then be written thanks to the linearity of Maxwell equations:

E ω n x; z X jmj<D∕λ Ẽi m E f m x; ze ik m x nd : (23)
This description working for cylindrical propagating beams, whose field distribution respects the 1D wave equation, can be generalized to 2D beams [START_REF] Chevalier | Electromagnetic modelization of spherical focusing on a one-dimensional grating thanks to a conical b-spline modal method[END_REF].

Construction of the SH Field

The following computation is similar to the previous construction of the fundamental field with some substantial differences: the first being the wavelength which is now equal to λ∕2. The first step consists of writing the complex second order polarization as a Fourier series:

P 2 D x; z X m 0 ∞ m 0 -∞ P2 m 0 ze ik m 0 x x ; (24)
where k m 0 x 2πm 0 ∕D and the wave vector is k k m 0 x ;k z such that jkj4π∕λ. All the modes, including the nonpropagative ones, are taken into account. The Fourier coefficients are given by

P2 m 0 z 1 D Z D∕2 -D∕2 P 2 D x; ze -ik m 0 x x dx: (25)
In Section 2.A.2, the nonlinear medium was decomposed into sub-layers generating fractions of the SH field from the source term P 2 x; z. Here each sub-layer generates a fraction of the SH field which is the sum of SH fields created from source terms for mode m 0 :P 2 m 0 ze -ik m 0 x x . The term dΠz now depends on m 0 :

dΠ m 0 z-iω P2 m 0 ze -ik m 0 x x dz: ( 26
)
The SH field is then computed in one period for every mode m 0 , thanks to the method described in Section 2.A.2. The latter is found in the other periods through the pseudo-periodicity condition written for this SH field: E shg m 0 x d; z E shg m 0 x; ze ik m 0 x d . The total field in the nth period eventually reads

E 2ω n x; z X m 0 ∞ m 0 -∞ E shg m 0 x; ze ik m 0 x nd : (27)
This result does not depend on the form of the nonlinear polarization, but is nonetheless limited to the undepleted pump regime. Propagation throughout the structure is also supposed to remain unperturbed by the presence of a nonlinear polarization term.

NUMERICAL RESULTS

A. Maker Fringes

Numerical simulations are performed in the case of Maker fringes in a 110 crystalline GaAs membrane (see Fig. 5), where a square sine function at a spatial frequency Δk k 2ω -2k ω 2n 2ω -n ω ω∕c is expected in the direction z of the incident plane wave [START_REF] Boyd | Nonlinear Optics[END_REF][START_REF] Maker | Effects of dispersion and focusing on the production of optical harmonics[END_REF]. The incident electric field is polarized along a 111 axis, and the following parameters are used: λ 2πc∕ω 5 μm, n 2 ω ε GaAs ω11.01, n 2 2ω ε GaAs 2ω 11.19 [START_REF] Kachare | Measurements of layer thicknesses and refractive indices in high-energy ion-implanted GaAs and GaP[END_REF] and χ 2 2∕ 3 p × 370 pm∕V (for polarization along the 111 axis) [START_REF] Boyd | Nonlinear Optics[END_REF]. The thickness of the crystal is equal to three times the coherence length L coh π∕Δk, and both faces have an anti-reflective coating at ω and 2ω, giving

jE 2ω zj 2 ∝ e -ik 2ω z sin Δkz 2 Δk - sin Σkz 2 Σk 2 ; ( 28
)
where Σk k 2ω 2k ω appears to be the frequency of the negligible oscillations which can be seen on the inset. They correspond to the sum of the contributions of upgoing waves (dU ).

A Gaussian focused beam at ω is now incident onto the same structure, with a waist W 5 μm and carrying a power of P 10 mW. The expected solution for an incoming 2D collimated Gaussian beam is a Bessel function multiplied by a Gaussian envelope [START_REF] Boyd | Nonlinear Optics[END_REF]. Here the result must be adapted to take 1D collimation into account, giving an approximate normalized SH field distribution:

E 2ω foc x; z ∝ Ax; z q 0 z p e -ik ω x 2 ∕q 0 z-i2k ω z ; (29)
where q 0 iπW 2 ∕λ, and

Ax; z iω ε 2 p c χ 2 Z z 0 e iΔkz 0 q 0 z 0 p dz 0 : (30)
Figure 6(c) shows the evolution of the square of the SH field absolute value at the center of the beam with respect to z. The waist of the beam, at the entrance of the medium, is equal to W 5 μm. The computation is performed over a depth equal to three times the coherence length of the SHG process. The analytical and numerical curves are found to be in good agreement. The difference is mainly because of the fact that the analytical solution relies on the hypothesis that the SH field has the same confocal parameter b k ω W 2 as the fundamental beam. The absolute values of the fundamental and SH fields are also displayed in Figs. 6(a) and 6(b) as functions of x and z to show the shape of the beam. These amplitudes are normalized to the maximum value of the field at ω.

B. SHG Exaltation in a Grating of Nanorods

As this method aims at designing nanostructures able to greatly enhance the nonlinear processes, a single grating of cylindrical nanorods is studied. The behavior of this architecture has already been investigated both in linear [START_REF] Gómez-Medina | Extraordinary optical reflection from sub-wavelength cylinder arrays[END_REF][START_REF] Ghenuche | Optical extinction in a single layer of nanorods[END_REF] and nonlinear cases [START_REF] Borisov | Second harmonic generation from arrays of sub-wavelength cylinders[END_REF]. These works highlighted a need for a strong linear exaltation of the electric field inside the cylinders. For a given wavelength λ and diameter R, such a nanostructure exhibits reflectivity resonances for peculiar values of the period d.At2πc∕ω 5 μm and R 0.2d, a resonance is found for a sub-wavelength period, d 4.36 μm. Other involved parameters are ε rods 2, θ 8.25°, χ 2 2 pm∕V incident field, and jE ω j1 × 10 6 V∕m. For the parameters Δz 0.01R and N 200 (corresponding to the situation Δz ∼ 0.002λ), the resulting field maps to be compared with Ref. [START_REF] Yuan | Analyzing second harmonic generation from arrays of cylinders using dirichlet-to-neumann maps[END_REF] are displayed in Figs. 7(a) and 7(b). Circles are fitted with crosses as drawn on the field maps: the cross-shaped section of the rods has total sides of length R and branches of length R∕4. The difference between these results and those of the Ref. [START_REF] Yuan | Analyzing second harmonic generation from arrays of cylinders using dirichlet-to-neumann maps[END_REF] can be explained by the fact that circular sections are identified with cross-shaped sections.

C. Convergence for a Grating of GaAs Nanorods

The convergence of the method is now investigated in the case of a grating of nanorods as a function of the number N of samples in a period along x and as a function of the discretization pace along z direction Δz. A grating (period d 4.9 μm) of GaAs nanorods with square section (R 0.2d) is studied. These nanorods are fabricated out of a 110 GaAs substrate and they are oriented along a 111 direction. The incoming light (2πc∕ω 5 μm) is TE-polarized and quasi-normally incident (θ 1°), and the amplitude of the incident electric field is jE ω j1 × 10 6 V∕m. First, to assess the convergence of the linear computed field, the relative error on the fundamental transmitted intensity is plotted in Fig. 8 as a function of N. Second, the relative error made on the value of the output zeroth order SH intensity as a function of N is displayed in Fig. 9 for Δz 10 -4 R, where I 2ω Δz 10 -4 R; N 1000 is taken as the reference value. This error decreases as N increases and goes under 10 -4 from N ≃ 200. This relative error is also displayed as a function of Δz for a fixed N 1000, with the same reference value. The same behavior is observed, highlighting the fact that the choices of Δz and N have to be balanced to obtain an accurate result.

CONCLUSION

A method for computing nonlinear effects in multilayered nanostructured media has been presented in this paper. On one hand, the method for plane waves works under the undepleted pump approximation and is efficient for a properly chosen z-sampling. It has been shown to converge with an increasing number of interpolations points and decreasing discretization pace along z. On the other hand, numerical simulations under a focused beam have been achieved by decomposing the incident light onto a plane wave basis. These results can be extended to 3D focused beams. Recent developments in the B-spline modal method have also introduced a Galerkin description which was shown to give even better convergence rates for linear problems [START_REF] Granet | Efficient implementation of b-spline modal method for lamellar gratings[END_REF].

Fig. 1 .Fig. 2 .

 12 Fig. 1. Typical setup of a grating layer (B) surrounded by two semiinfinite media (A, C) separated by two interfaces α and β. It is illuminated by a TE-polarized plane wave E. The incident wave vector k makes an angle θ with the optical axis z in the xOz plane. The structure is infinite along the y direction.

Fig. 3 .

 3 Fig.3. Scheme of the reflections and propagations undergone by the upgoing waves dU (resp. the downgoing waves dU -) in the nonlinear medium from height z (resp. z -). The sub-layer of thickness dz at height z generates a source term dΠz-iωP 2 zdz.

Fig. 4 .

 4 Fig. 4. Scheme of a cylindrical collimated beam. It is infinite along the y axis.

Fig. 5 .

 5 Fig. 5. Absolute square value of the SH electric field, normalized by the maximum value of the analytical curve for Maker fringes in 110 crystalline GaAs. The inset zooms in the oscillations at frequency Σk. Involved parameters are 2πc∕ω 5 μm, ε GaAs ω11.01, and ε GaAs 2ω11.19, θ 0°.

Fig. 6 .

 6 Fig. 6. Fields maps (a) of the focused fundamental field and (b) of the SH field in 1 10 crystalline GaAs, normalized to the maximum amplitude of the fundamental field. (c) Absolute square value of the SH electric field at x 150 μm as a function of z, normalized to the maximum value of the analytical curve. The waist W of the Gaussian beam is located at the entrance of the crystal and is equal to 5 μm. Involved parameters are 2πc∕ω 5 μm, ε GaAs ω11.01, ε GaAs 2ω11.19, χ 2 2∕ 3 p × 370 pm∕V, and θ 0°. The beam carries a power P 10 mW.

Fig. 7 .

 7 Fig. 7. Fundamental and SH absolute values of the electric field in the plane perpendicular to the rod axis of a grating of GaAs nanorods. The electric field is polarized along the rods. Involved parameters are period d 4.36 μm, R 0.2d, 2πc∕ω 5 μm, χ 2 2 pm∕V, ε rods 2, θ 8.25°, and incident field jE ω j1 × 10 6 V∕m. The cross-shaped section of the rods has total sides of length R and branches of length R∕4.

Fig. 8 .

 8 Fig. 8. Relative error on the value of the output fundamental intensity as a function of the number of interpolation points N, with the reference value I ω N 100017.33 GW∕m 2 . The structure is a grating of GaAs nanorods; involved parameters are period d 4.90 μm, R 0.2d, 2πc∕ω 5 μm, ε GaAs ω11.01, ε GaAs 2ω11.19, θ 1°, and incident field jE ω j1 × 10 6 V∕m.

Fig. 9 .

 9 Fig. 9. Relative error on the value of the output zeroth order SH intensity in the zero diffractive order (a) as a function of the number N of samples along x in a period for Δz 10 -4 R, and (b) as a function of Δz, the discretization pace along z direction for N 1000. The reference value is I 2ω Δz 10 -4 R; N 10001.817757 W∕m 2 . The structure is a grating of GaAs nanorods; involved parameters are the same as those in Fig. 8 with χ 2 150 pm∕V.


Tarek Ahmed-Ali
email: tarek.ahmed-ali@ensicaen.fr

Emilia Fridman
email: emilia@eng.tau.ac.il

Fouad Giri
email: fouad.giri@unicaen.fr

Laurent Burlion
email: laurent.burlion@onera.fr

Françoise Lamnabhi-Lagarrigue

& Lamnabhi-Lagarrigue

Using exponential time-varying gains for sampled-data stabilization and estimation

This paper provides exponential stability results for two system classes. The first class includes a family of nonlinear ODE systems while the second consists of semi-linear parabolic PDEs. A common feature of both classes is that the systems they include involve sampled-data states and a time-varying gain. Sufficient conditions ensuring global exponential stability are established in terms of Linear Matrix Inequalities (LMIs) derived on the basis of Lyapunov-Krasovskii functionals. The established stability results prove to be useful in designing exponentially convergent observers based on sampled-data measurements. It is shown throughout simulated examples from the literature that the introduction of time-varying gains is beneficial to the enlargement of sampling intervals while preserving the stability of the system.

Introduction

Designing sampled-data observers and controllers has been a hot topic in recent years, see e.g. [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF] and reference list therein. In this regard, a long standing issue is how to enlarge the sampling intervals [START_REF] Heemels | An introduction to eventtriggered and self-triggered control[END_REF] while ensuring global exponential stability. In a recent paper [START_REF] Cacace | Nonlinear systems with multiple timevarying measurement delays[END_REF], it has been shown that the introduction of time-varying gains in a specific class of observers improves their exponential convergence properties in presence of measurement delay. Presently, these properties are investigated in presence of measurement sampling. To this end, we consider two classes of sampled-data systems and analyze their exponential stability. The considered classes are respectively consisting of nonlinear globally Lipschitz ODEs and semi-linear parabolic PDEs. A common feature of both classes is that the systems they include are allowed to involve a time-varying gain of the form e -η(t-t k ) with η > 0 a tuning parameter, where t k (k = 0, 1, . . .) are sampling instants. It turns out that, the first family, including ODE systems, is a generalization of that dealt with in [START_REF] Cacace | Nonlinear systems with multiple timevarying measurement delays[END_REF]. For both classes of systems, we establish sufficient conditions for global exponential stability in terms of Linear Matrix Inequalities (LMIs) derived from Lyapunov-Krasovskii functionals. Then, it is shown that these stability results are useful in designing sampled-data observers with timevarying gains. As the established LMIs conditions involve both the tuning parameter η and the maximum sampling interval h, these parameters can then be used to improve the observer convergence properties. Actually, it is checked through several simulated examples that the utilization of the above time-varying gain entails significant enlargement of the maximum sampling interval, compared with the constant gain case (corresponding to η = 0). It is worth noting that, the present theoretical stability results can also be used in sampled-data control design improving exponential stability properties and enlarging sampling intervals. A part of the present results, namely those concerning ODEs, have been presented in our conference paper [START_REF] Ahmed-Ali | A new approach to enlarging sampling intervals for some sampled-data systems and observers[END_REF].

The paper is organized as follows: in Section 2, the first stability result, concerning nonlinear ODE systems, is stated and applied to observer design; in Section 3, the second stability result, concerning semi-linear PDE systems, is stated and applied to observer design; a conclusion and reference list end the paper. Some technical proofs are appended.

Notations and preliminaries

Throughout the paper the superscript T stands for matrix transposition, R n denotes the n-dimensional Euclidean space with vector norm |.|, R n×m is the set of all n × m real matrices, and the notation P > 0, for P ∈ R n×n , means that P is symmetric and positive definite. In Symmetric matrices, symmetric terms are denoted * ; λ min (P) (resp. λ max (P)) denotes the smallest (resp. largest) eigenvalue. The notation (t k ) k≥0 refers to a strictly increasing sequence such that t 0 = 0 and lim k→∞ t k = ∞. The sampling periods are bounded i.e. 0 < t k+1t k < h for some scalar 0 < h < ∞ and all k = 0, 1, . . . , ∞. We also define the variable τ (t

) = t -t k , t ∈ [t k , t k+1

Sampled-data globally Lipschitz nonlinear ODEs

System description and stability result

We are considering a class of sampled-data nonlinear systems described by the following equation:

ẋ(t) = A 0 x(t) + A 1 e -η(t-t k ) x(t k ) + φ(x(t)), t ∈ [t k , t k+1 ) (1)
where x(t) ∈ R n ; the scalar η ≥ 0; A 0 , A 1 are constant matrices with appropriate dimensions. As in Bar Am and [START_REF] Fridman | Systems and control: foundations and applications, Introduction to time-delay systems: analysis and control[END_REF], the function φ is supposed to be class C 1 with uniformly bounded Jacobian φ x , satisfying φ(0) = 0 and

φ T x (x)φ x (x) ≤ M ∀x (2)
for some positive constant n × n-matrix M. Using Jensen's inequality it is readily checked that (2) implies the following inequality:

 1 0 φ T x (sx)ds  1 0 φ x (sx)ds ≤ M.
Remark 1. Eq. ( 1) may represent a networked control system described by

ẋ(t) = A 0 x(t) + φ(x(t)) + Bu(t),
with the communication network placed between the sensor and the controller (but no network is placed between the controller and the actuator). Assuming that the discrete-time state measurements

x(t k ) are transmitted through the communication network from the sensor to controller, consider the state-feedback,

u(t) = e -η(t-t k ) Kx(t k ), t ∈ [t k , t k+1 ),
where K is a gain and η > 0 is a scalar. It turns out that the resulting closed-loop system fits Eq. ( 1) with A 1 = BK .

As in [START_REF] Cacace | Nonlinear systems with multiple timevarying measurement delays[END_REF], introduce the following change of coordinates z(t) = e ηt x(t) with η > 0. Then one gets

ż(t) = ηz(t) + A 0 z(t) + A 1 z(t k ) +  1 0 φ x (sx(t))ds  e ηt x(t), t ∈ [t k , t k+1 ) (3)
which is rewritten as follows:

ż(t) = (ηI n + A 0 ) z(t) + A 1 z(t k ) +  1 0 φ x (sx(t))ds  z(t), t ∈ [t k , t k+1 ). (4)
Following [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF], consider the following Lyapunov-Krasovskii functional for (4):

V (t) = V (t) + V X (t) (5)
with

V (t) = z T (t)Pz(t) + (t k+1 -t)  t t k żT (s)U ż(s)ds, P > 0, U > 0, t ∈ [t k , t k+1 ) and V X (t) = (t k+1 -t)ξ T     X + X T 2 -X + X 1 * -X 1 -X T 1 + X + X T 2     ξ ,
where ξ (t) = col{z(t), z(t k )}, X and X 1 are n × n matrices. The positiveness of ( 5) is ensured if the following LMI holds [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF]:

    P + h X + X T 2 hX 1 -hX * -hX 1 -hX T 1 + h X + X T 2     > 0. (6)
Using the definition of z(t), we can see that the exponential stability of system (1) is guaranteed if:

V (t) + 2αV (t) ≤ 0 t ∈ [t k , t k+1 ) (7)
for some scalar α ∈ (-η, 0] (note that the scalar α is allowed to be negative). Indeed, if (7) is satisfied one has,

V (t) ≤ -2αV (t) =⇒ |z(t)| ≤   V |t=0 √ λ min (P)  e -αt .
Then, using the fact that z(t) = e ηt x(t), one gets:

|x(t)| ≤   V |t=0 √ λ min (P)  e -(η+α)t .
From the above inequality, one sees that the exponential convergence is guaranteed if η + α > 0. Since the parameter η is positive and free, it is sufficient to let α ∈ (-η, 0] for ensuring an exponential convergence with a decay rate η + α. In the following proposition, it is shown that the property (7), and resulting exponential stability with a decay rate η + α > 0, are actually ensured under well established sufficient conditions, expressed in terms of LMIs.

Proposition 1. Consider the system (1) with possibly varying sampling-intervals subject to t k+1t k ≤ h with some scalar h > 0.

Given η > 0 and α ∈ (-η, 0], let there exist n × n matrices P > 0, U > 0, X , X 1 , P 2 , P 3 , T , Y 1 , Y 2 and a scalar λ > 0 that satisfy the LMI (6) and the following LMIs:

Ψ (0)     Φ 11 -X α Φ 12 + X τ (t) Φ 13 + X 1α P T 2 * Φ 22 + hU Φ 23 -X 1τ (t) P T 3 * * Φ 33 -X 2α 0 * * * -λI n       τ(t)=0 < 0 (8) and Ψ (h)       Φ 11 -X α| Φ 12 Φ 13 + X 1α hY T 1 P T 2 * Φ 22 Φ 23 hY T 2 p T 3 * * Φ 33 -X 2α hT T 0 * * * -hU 0 * * * * -λI n         τ(t)=h < 0, ( 9
)
where

Φ 11 = (A 0 + ηI n ) T P 2 + P T 2 (A 0 + ηI n ) + 2αP -Y 1 -Y T 1 + λM Φ 12 = P -P T 2 + (A 0 + ηI n ) T P 3 -Y 2 Φ 13 = Y T 1 + P T 2 A 1 -T Φ 22 = -P 3 -P T 3 Φ 23 = Y T 2 + P T 3 A 1 Φ 33 = T + T T X τ (t) = (h -τ (t)) X + X T 2 X α = (1 -2α(h -τ (t))) X + X T 2 X 1τ (t) = (h -τ (t))(X -X 1 ) X 1α = (1 -2α(h -τ (t)))(X -X 1 ) X 2α = (1 -2α(h -τ (t))) X + X T -2X 1 -2X T 1 2
.

Then, the system (1) is exponentially stable with a decay rate η+α. Moreover, if the above LMIs hold with α = -η, then the system is exponentially stable with a small enough decay rate in (0, ϵ) for some ϵ ∈ (0, -α).

See Appendix A for the proof.

Remark 2. Practically, the following procedure is used to get suitable values of the design parameters involved in Proposition 1.

The search procedure is started by taking a small value of η (e.g. η = 0). Then, one proceeds as follows:

(1) Let α = -η.

(2) Take a small value of h and check the LMIs.

(3) Increase h until the LMIs are no longer feasible. Then, retain the value of h.

(4) Increase the value of η and repeat Steps (1)-( 3).

(5) The above procedure is stopped when the increase of η does not entail an increase of h.

Remark 3. Proposition 1 clearly shows that the maximum sampling interval h that preserves the exponential stability is made dependent on the design parameters η and α. Additional highlights will be provided in Example 1 (given hereafter) where it is shown that the maximum h preserving stability may be significantly enlarged by tuning η. On the other hand, Proposition 1 also con- firms the somewhat obvious fact that the exponential decay rate η + α < η is also depending on η and α.

Remark 4. In the linear case, the function φ (see ( 1)) is zero and M = 0 in (2). Then, it is sufficient to verify the feasibility of the smaller-order λ-independent LMIs ( 8) and ( 9) with the deleted last column and row. Indeed, the feasibility of the smaller-order LMIs implies the feasibility of the full-order

ẋ(t) = A 0 x(t) + A 1 e -η(t-t k ) x(t k ), t ∈ [t k , t k+1 )
LMIs with some λ > 0.

Remark 5. If in (1) the matrix A 1 = BK depends on the unknown gain K , the inequalities ( 8) and ( 9) become nonlinear. However, by using the standard arguments for the controller design via the descriptor method [START_REF] Fridman | Systems and control: foundations and applications, Introduction to time-delay systems: analysis and control[END_REF], where it is assumed that P 2 = εP 3 with a tuning parameter ε, one can easily derive LMIs for finding K .

Example 1

To illustrate the result of Proposition 1, consider the following system

ẋ(t) = -x(t) + sin(x(t)) -Ke -η(t-t k ) x(t k ), t ∈ [t k , t k+1 ) (10)
with K = 1. Here we have A 0 = -1, A 1 = -K = -1, and M = 1. Then, application of Proposition 1 with α = -η with various values of η leads to different values of h that preserve the exponential stability of the system (see Table 1). Thus, (10) with η = 1 achieves exponential stability for more than twice larger sampling interval compared with the constant gain (i.e. case η = 0). As a matter of fact, the parameter η cannot be infinitely increased. Presently, the maximum value of η that still yields an increase of the sampling interval is η = 1.

Example 2

Consider the following much studied system (see e.g. [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF][START_REF] Zhang | Stability of networked control systems[END_REF]: -[3.75, 11.5]. It is well-known that the above system with η = 0 remains stable under constant sampling with h < 1.72 and becomes unstable for constant sampling periods with h > 1.73. Applying the result of Proposition 1, one obtains for η = 0.8 and α = -0.8 a maximum sampling interval value of h = 2.43 that preserves the stability.

ẋ(t) =  0 1 0 -0.1  x(t) + e -η(t-t k )  0 0.1  Kx(t k ) t ∈ [t k , t k+1 ) with K =

Example 3

Consider the following closed loop system:

ẋ(t) =  0 1 -1 -2  x(t) +  0 1  e -η(t-t k ) Kx(t k ) t ∈ [t k , t k+1 ) (11) with Ke -η(t-t k ) =  -1 1  e -η(t-t k )
is the regulator gain. Applying the procedure of Remark 2 to system (11), one obtains the results of Table 2. This shows that the time-varying nature of the regulator gain is actually beneficial to the enlargement of the sampling interval. The analysis of the closed loop system (11) in the case of a constant regulator gain (i.e. η = 0) has been made in several previous studies using various approaches yielding different maximum sampling intervals. A comparison of the obtained maximum sampling intervals is made in [START_REF] Kao | On robust stability of aperiodic sampled-data systemsan integral quadratic constraint approach[END_REF], see Example 2 in Section 5 therein. The highest maximum sampling interval, obtained in [START_REF] Kao | On robust stability of aperiodic sampled-data systemsan integral quadratic constraint approach[END_REF], is h = 1.728.

Clearly, this value is much smaller than the maximum sampling interval achieved by using a time-varying regulator gain. The above examples illustrate the fact that, when using a time varying gain, the maximum sampling interval h that preserves the stability can be much larger compared to the constant gain case (η = 0).

Application to sampled-data observer design

Consider the class of nonlinear systems

 ẋ(t) = Ax(t) + f (x(t)) y(t) = Cx(t) (12)
where x ∈ R n and A, C are matrices with appropriate dimensions.

It is supposed that the function f satisfies assumption (2) and y is accessible to measurements only at instants t k . We also suppose that the pair (A, C ) is detectable. The following observer is then proposed:

ẋ(t) = Ax(t) + f (x(t)) -Ke -η(t-t k ) (C x(t k ) -y(t k )) t ∈ [t k , t k+1 ), (13)
with η > 0 and K is a matrix gain of appropriate dimension. Our goal is to determine a matrix K and a maximum sampling period h so that, the observation error x(t) = x(t)x(t) is globally exponentially stable. It is readily checked that this error satisfies the following equation:

ẋ(t) = Ax(t) + f (x(t)) -f (x(t)) -Ke -η(t-t k ) (C x(t k ) -y(t k )) t ∈ [t k , t k+1 ).
Considering the change of coordinates z 1 (t) = e ηt x(t), one gets the following equation

ż1 (t) = (ηI n + A)z 1 (t) -KCz 1 (t k ) + e ηt (f (x(t)) -f (x(t))). ( 14
)
Using the fact that,

f (x(t)) -f (x(t)) =  1 0 f x (x(t) + s(x(t) -x(t)))ds  x(t),
Eq. ( 14) can be rewritten as follows:

ż1 (t) = (ηI n + A)z 1 (t) + A 1 z 1 (t k ) +  1 0 f x (x(t) + s(x(t) -x(t)))ds  z 1 (t), t ∈ [t k , t k+1 ) (15)
with A 1 = -KC. Clearly, Eq. ( 15) falls in the class of systems described by Eq. ( 4). Therefore, the result of Proposition 1 can directly be applied to get sufficient conditions for the observer (13) to be exponentially convergent. This is illustrated in the next example.

Example 4

Consider the following nonlinear system

 ẋ1 (t) = x 2 (t) ẋ2 (t) = -2x 1 (t) -x 2 (t) + 0.2 sin(x 2 (t)) y(t) = Cx(t) = x 1 (t).
Clearly, this system is of the form (12 with M = diag{0, 0.04}. Then the observer ( 13) is written as follows:

) with A =  0 1 -2 -1  , C = (1, 0), and f (x(t)) =  0 0.2 sin(x 2 (t))  T .
   ẋ1 (t) = x2 (t) -K 1 e -η(t-t k ) (x 1 (t k ) -y(t k )) ẋ2 (t) = -2x 1 (t) -x2 (t) + 0.2 sin(x 2 (t)) -K 2 e -η(t-t k ) (x 1 (t k ) -y(t k ))
where

K = [K 1 K 2 ] and A 1 = -KC =  -K 1 0 -K 2 0



. Applying Proposition 1 with K = [1.65, -0.85] T and α = -η, one obtains the results of Table 3. This clearly shows that the maximum sampling interval h can be significantly increased thanks to the time-varying gain feature of the observer. Remark 6. Note that high gain observers are particular cases of the class of observers described by (13). They correspond to K = θ ∆ -1 L where θ ≥ 1, ∆ = diag{1, . . . , 1/θ n-1 } and L is such that A -LC is Hurwitz. Then the function in ( 12) must be triangular [START_REF] Gauthier | A simple observer for nonlinear systems applications to bioreactors[END_REF]. Compared to the timevarying gain observer proposed in [START_REF] Farza | Continuous-discrete time observers for a class of mimo nonlinear systems[END_REF], our observer (13) enjoys a simpler time-varying gain (because it is scalar) and offers more flexibility because the parameter η is independent on the gain K . In this respect, one might also note that the timeinvariant-gain version of the observer proposed in [START_REF] Farza | Continuous-discrete time observers for a class of mimo nonlinear systems[END_REF] is nothing other than an obvious rewriting in impulsive form of the original observer of [START_REF] Karafyllis | From continuous-time design to sampleddata design of nonlinear observers[END_REF].

Sampled-data semi-linear parabolic PDEs

System description and stability result

In this section, it is shown that results similar to those of Section 2 can be established for infinite-dimensional systems. Let us consider the class of systems governed by the following partial differential equation:

u t (x, t) = u xx (x, t) + φ(x, u(x, t), t)u(x, t) -Ke -η(t-t k ) u(x j , t k ) t ∈ [t k , t k+1 ), x ∈ [x j , x j+1 ), (16) with xj =
x j+1 +x j 2 (j = 0, . . . , N -1), where the points x j divide the interval [0, l]

such that 0 = x 0 < • • • < x N = l, under Neumann boundary conditions u x (l, t) = u x (0, t) = 0, t ≥ 0. ( 17
)
We also suppose that the function φ is C 1 and is possibly unknown but φ m ≤ φ ≤ φ M , where φ m and φ M are known bounds. The sampling interval in space may be variable but bounded

x j+1 -x j ≤ ∆ ≤ l.
By arguments developed in [START_REF] Fridman | Robust sampled-data control of a class of semilinear parabolic systems[END_REF], there exists a unique strong solution of ( 16) initialized by

u(•, 0) ∈ H 1 (0, l) : u x (0, 0) = u x (l, 0) = 0.
η > 0, h > 0, R > 0, 2η > δ 1 > 0, 2δ ∈ (δ 1 -2η, 0]. Choose K > φ M + η -π 2
l 2 and suppose that there exist scalars p i > 0, r > 0, and y i (i = 1, 2, 3), satisfying the inequality

KR -1 (p 2 + p 3 ) ≤ π δ 1 p 3 (19)
and the following LMIs:

Φ i |φ=φ m < 0, Φ i |φ=φ M < 0, i = 0, 1 (20
χ 11 = 2δp 1 + 2p 2  φ + η + KR 2π  -2y 1 χ 12 = p 1 -p 2 + p 3 (φ + η) -y 2 χ 13 = y 1 -y 3 -Kp 2 χ 22 = -2p 3 + KRp 3 π χ 23 = y 2 -Kp 3 λ = 2 π 2 l 2 (p 2 -δp 3 ) χ 33 = 2y 3 -δ 1 p 1 .
Then, the unique strong solution of the system (16) initiated by (18) satisfies the following bound, for all t ≥ 0:

 l 0 [p 1 u 2 (x, t) + p 3 u 2 x (x, t)]dx ≤ e (-2(δ+η)+δ 1 )t  l 0 [p 1 u 2 (x, 0) + p 3 u 2 x (x, 0)]dx, ( 21
)
accordingly the system (16) is exponentially stable with a decay rate (δ + η)δ 1 /2. Moreover, if the LMIs (20) hold with 2δ = δ 1 -2η, then the system (16) is exponentially stable with a decay rate in (0, ϵ) for some ϵ ∈ (0, δ 1 ).

See the proof placed in Appendix B.

Remark 7. A practical selection procedure similar to the ODE case (Remark 2) is used to find out suitable values of the design parameters of Proposition 2. Again, the search is started by considering a small value of η, (e.g. η = 0).

(1) Choose δ 1 arbitrarily in the interval (0, 2η) and let δ = δ 1 /2-η.

(2) Take a small value of h and check the LMIs.

(3) Increase h until the LMIs are no longer feasible. Then, retain the value of h.

(4) Increase the value of η and repeat Steps (1)-( 3).

(5) The above procedure is stopped when the increase of η entails no increase of h.

Example 5 Consider the following system:

u t (x, t) = u xx (x, t) -Ke -η(t-t k ) u(x j , t k ) x ∈ [x j , x j+1 ), t ∈ [t k , t k+1 )
under the Neumann boundary conditions. This system fits system structure ( 16), with φ(x, u) = 0. Applying Proposition 2 with l = π , K = 1, ∆ = π /7, R = 1, one gets for different values of η the maximum sampling interval that preserves the exponential stability. The results thus obtained are summarized in Table 4 which shows that the maximum sampling interval h can be substantially enlarged by tuning η.

Application to sampled-data observer design

Consider the following semi-linear diffusion equation:

u t (x, t) = u xx (x, t) + f (u(x, t), x, t) (22)
with Neumann condition u x (0, 0) = u x (l, 0) = 0. The system output, y(t k ) = u(x j , t k ), for some j, is only available at sampling

instants t k . It is supposed that the function f is known, of class C 1 ,
and satisfying f m ≤ f u ≤ f M , for some scalar constants f m and f M .

The following observer structure is considered:

ût (x, t) = ûxx (x, t) + f (û(x, t), x, t) -Le -η(t-t k ) (û(x j , t k ) -y(t k )) x ∈ [x j , xj+1 ), t ∈ [t k , t k+1 ) (23)
with ûx (0, t) = ûx (l, t) = 0. It is readily checked that the observation error e(x, t) = û(x, t)u(x, t) undergoes the following equation:

e t (x, t) = e xx (x, t) + φ(e(x, t), x, t)e(x, t) -Le -η(t-t k ) e(x j , t k ), x ∈ [x j , xj+1 ), t ∈ [t k , t k+1 ) (24)
where φ =  l 0 f u (û + θ e, x, t)dθ , with boundary conditions e x (0, t) = e x (l, t) = 0. Clearly, Eq. ( 24) is of the form ( 16). There- fore, one can make use of Proposition 2 to determine the observer gain L as well as a suitable value of the maximum h, while guaranteeing the exponential convergence of the observer (23).

Example 6

Heat/mass transfer systems with heat generation or volumetric chemical reactions, e.g. chemical tubular reactor, are examples of systems that can be described by ( 16), see [START_REF] Boskovic | Backstepping control of chemical tubular reactor[END_REF] and references therein. Indeed, such systems undergo a heat equation with non-constant coefficient i.e.

u t = ϵ u u xx + λ αβ (x)u(x, t), x ∈ (0, 1) (25) with u x (0, t) = u x (1, t) = 0, where λ αβ (x) = 2 cosh 2 (αx -β) . (26)
Let ϵ u = 1, α = 4 and β = 2 and consider the following observer:

ût = ûxx + λ αβ (x)û(x, t) -Le -η(t-t k ) (û(x j , t k ) -y(t k )), x ∈ (0, 1) (27)
with ûx (0, t) = ûx (1, t) = 0. It is readily checked that the observation error e(x, t) = û(x, t)u(x, t) undergoes Eq. ( 24) with φ(e, x, t) = φ(x) = λ αβ (x) and the boundary conditions e x (0, t) = e x (1, t) = 0. It can be shown that one presently has φ M = 30 and φ m = 0. Applying the above procedure with R = 1, ∆ = 0.001, L = 0.11 + φ M + ηπ 2 , several possible sampling interval values, corresponding to different values of the parameter η, have been obtained, see Table 5. It is readily seen that the observer time-varying gain associated with η = 2 provides a maximum sampling interval that is nearly 10 times larger than the case of fixed observer gain (i.e. case η = 0).

Conclusion

Novel stability results, stated in Propositions 1 and 2, are established for two classes of systems respectively described by ODEs and PDEs. Both studied classes involve the exponentially decaying term e -η(t-t k ) and all stability conditions are expressed in terms of LMIs. The new results conditions are shown to be useful in designing sampled-data observers with exponentially decaying gains. Furthermore, it is checked through several numerical examples that the introduction of exponentially decaying gains may lead to a substantial enlargement of the maximum sampling intervals while preserving stability. A theoretical proof that, timevarying gain entails sampling interval enlargement is yet to be found. Finally, the extension of the results presented in the present paper, to the case of disturbances is a topic for our future research.

Appendix A. Proof of Proposition 1

Let α ∈ (-η, 0] and differentiate V (t). After some simple computations one obtains:

V (t) + 2αV (t) + λz T (t) ×  M -  1 0 φ T x (sx)ds  1 0 φ x (sx)ds  z(t) ≤ 2z T (t)P ż(t) + (t k+1 -t)ż T (t)U ż(t) -  t t k żT (s)U ż(s)ds -ξ T     X + X T 2 -X + X 1 * -X 1 -X T 1 + X + X T 2     ξ + (t k+1 -t)[ż T (t)(X + X T )z(t) + 2ż T (t)(-X + X 1 )z(t k )] + 2αz T (t)Pz(t) + 2α(t k+1 -t)ξ T     X + X T 2 -X + X 1 * -X 1 -X T 1 + X + X T 2     ξ + λz T (t)Mz(t) -λz T (t)  1 0 φ T x (sx)ds  1 0 φ x (sx)ds  z(t) (A.1)
using the well known Jensen's inequality, one gets,

 t t k żT (s)U ż(s)ds ≥ (t -t k )v T 1 Uv 1 with v 1 = 1 t-t k  t t k
żT (s). According to the descriptor approach [START_REF] Fridman | New Lyapunov Krasovskii functionals for stability of linear retarded and neutral type systems[END_REF], the left-hand sides of the equations,

0 = 2  z T (t)Y T 1 + żT (t)Y T 2 + z T (t k )T T  × [-z(t) + z(t k ) + (t -t k )v 1 ] 0 = 2  z T (t)P T 2 + żT (t)P T 3  ×  (A 0 + ηI n )z(t) + A 1 z(t k ) +  1 0 φ x (sx)ds  z(t) - ż (A.2)
can be added to V (t) + 2αV (t) where Y 1 , Y 2 , T , P 2 , P 3 are free matrices. Let us define the augmented vector,

µ = col  z(t), ż(t), z(t k ), v 1 ,  1 0 φ x (sx)ds  z(t)
 .

Then, combining (A.2) and (A.1), one gets:

V (t) + 2αV (t) + λz T (t)  M -  1 0 φ T x (sx)ds  1 0 φ x (sx)ds  z(t) ≤ µ T Ψ µ < 0 (A.3) provided that Ψ       Φ 11 -X α Φ 12 + X τ (t) Φ 13 + X 1α τ (t)Y T 1 P T 2 * Φ 22 + s.U Φ 23 -X 1τ (t) τ (t)Y T 2 p T 3 * * Φ 33 -X 2α τ (t)T T 0 * * * -τ (t)U 0 * * * * -λI n       < 0. (A.4)
with s = t k+1t. Writing the last matrix inequality for τ (t) → 0 and τ (t) → h leads to LMIs ( 8) and ( 9), respectively. Now, using arguments of Lemma 2 in [START_REF] Fridman | A refined input delay approach to sampled-data control[END_REF], it is not difficult to see that, if the LMIs ( 6), ( 8) and ( 9) are feasible for some h > 0, then they are also feasible for all h ∈ (0, h]. Now, consider the case where α = -η. If the LMIs (6), ( 8) and ( 9) are feasible, then by continuity, they are also feasible with the same η and ᾱ = -η+ϵ for small enough ϵ > 0. Then, an inequality like (A.3) holds with α being replaced by ᾱ. Then, one gets

V (t) + ᾱV ≤ 0 which yields |x(t)| ≤   V |t=0 √ λ min (P)  e -(η+ ᾱ)t =   V |t=0 √ λ min (P)  e -ϵt .
This completes the proof of Proposition 1.

Appendix B. Sketch of proof of Proposition 2

Let us introduce the change of coordinates ξ (x, t) = e ηt u(x, t) with η > 0. Then it is checked that ξ (x, t) undergoes the equation:

ξ t (x, t) = ξ xx (x, t) + (φ(x, u(x, t), t) + η)ξ (x, t) -K ξ (x j , t k ) t ∈ [t k , t k+1 ), x ∈ [x j , xj+1 ). (B.1)
Consider the following Lyapunov-functional inspired by Fridman and Blighovsky (2012):

V 1 (t) = p 1  l 0 ξ (x, t) 2 dx +  l 0  p 3 ξ 2 x (x, t) + r(t k+1 -t) ×  t t k e 2δ(s-t) z 2 s (x, s)ds  dx, t ∈ [t k , t k+1 ).
Notice that the above functional is continuous i.e. it satisfies

V 1 (t k ) = V 1 (t - k ).
V1 (t) + 2δV 1 (t) -δ 1 V 1 (t k ) ≤  l 0 η T 2 Φs η 2 dx +  ∆ π KR -1 (p 2 + p 3 ) -δ 1 p 3   l 0 ξ 2 x (
V1 (t) + 2δV 1 (t) -δ 1 V 1 (t k ) ≤ 0 t ∈ [t k , t k+1 ).
Since 2δ < 0, one easily gets | V1 (t)| ≤ (-2δ + δ 1 ) sup

t k ≤s<t (V 1 (s)) t ∈ [t k , t k+1 ).
Now, let us define the following continuous functional: Then, by using the comparison Lemma 2.12 in [START_REF] Karafyllis | Stability and stabilization of nonlinear systems[END_REF], one obtains 20) are feasible, then by continuity, they are also feasible with the same η and δ1 = 2δ + 2ηϵ for small enough ϵ > 0. Then, an inequality like (B.4) holds with δ 1 being replaced by δ1 . Then, one gets

ψ(t) = sup t k ≤s<t (V 1 (s)), t ∈ [t k , t k+1
V 1 (t) ≤ sup t k ≤s<t (V 1 (s)) ≤ V 1 (0)e (-2δ+δ
V 1 (t) ≤ sup t k ≤s<t (V 1 (s)) ≤ V 1 (0)e (-2δ+ δ1 )t , t ≥ 0 (B.5) which yields  l 0 [p 1 u 2 (x, t) + p 3 u 2 x (x, t)]dx ≤ V 1 (0)e -ϵt .
This completes the proof of Proposition 2.

Proposition 2 .

 2 Consider the class of systems described by (16) and the positive scalars

  ). H 1 (0, l) is the Sobolev space of absolutely continuous functions z : (0, l) → R with the square integrable derivative d

dx . Given a two-argument function u(x, t), its partial derivatives are denoted u t = ∂u ∂t , u xx = ∂ 2 u ∂x 2 .

Table 1

 1 Max. values of h vs. η preserving the stability of Example 1.

	η	0	0.5	0.7	1
	h	1.38	1.84	2.36	3.17

Table 2

 2 Max. values of h vs. η preserving the stability of Example 3.

	η	0.2	0.6	1	1.5
	h	2.03	3.38	4.69	5.60

Table 3

 3 Max. values of h vs. η preserving the stability of Example 4.

	η	0	1	2	2.4
	h	0.99	1.64	2.78	2.86
		It is readily checked that the function f satisfies inequality (2)

  )

	where
	Φ 0 Φ 1	χ 12 hr + χ 22 χ 23 χ 13 χ 33 * χ 11 -λ χ 12 χ 13 hy 1   χ 11 -λ * *    χ 22 χ 23 hy 2 * * * * * -hr * χ 33 hy 3	  
	with

Table 4

 4 Max.values of h vs. η that preserve the stability of Example 5.

	η	0	0.25	0.75	1	2
	δ 1	0.5	0.5	0.5	0.5	0.5
	δ h	0.25 1.27	0 1.95	-0.25 1.8	-0.5 1.73	-0.75 1.41

Table 5

 5 Max. values of h vs. η that preserve the stability of Example 6.

	η	0	0.5	0.75	1	2
	δ 1	0.5	0.5	0.5	0.5	0.5
	δ h × 10 3	0.25 3	-0.25 16.5	-0.5 19.9	-0.75 22.6	-1.75 29

  The rest of the proof matches mutatis mutandis a similar proof in[START_REF] Fridman | Robust sampled-data control of a class of semilinear parabolic systems[END_REF]. Relying on the descriptor technique, the two following equalities are considered:[p 2 ξ (x, t) + p 3 ξ t (x, t)] × [-ξ t (x, t) + ξ xx (x, t) + (φ(x, u(x, t)) + η)ξ (x, t) -K ξ (x, t k )]dx [p 2 ξ (x, t) + p 3 ξ t (x, t)] [y 1 ξ (x, t) + y 2 ξ t (x, t) + y 3 ξ (x, t k )] × [-ξ (x, t) + ξ (x, t k ) + (tt k )v 2 (x, t)] dx (x, s)dswhere the first equality is directly got from (B.1) while the second is obtained using Leibnitz-Newton formula. Adding the right sides of equalities (B.2), (B.3) to V1 (t) + 2δV 1 (t)δ 1 V 1 (t k ), one obtains after some computations the following inequality:

			 l
	0 = 2	0
		+ 2K	N-1  j=1	 x j+1 x j
			 x
		×  l xj	ξ ζ (ζ , t k )dζ dx	(B.2)
	0 = 2	0
					(B.3)
	with
	v 2 =	1 t -t k	 t t k	ξ s

  x, t k )dx with η 2 = col{ξ (x, t), ξ t (x, t), ξ (x, t k ), v 2 } t k )y 1 * (t k+1t)r + χ 22 χ 23 (tt k )y 2Just as in[START_REF] Fridman | Robust sampled-data control of a class of semilinear parabolic systems[END_REF], it can be seen that if the LMIs (20) are feasible, then Φs ≤ 0. If further (19) is fulfilled, then the following inequality holds:

	and
	Φs =	   χ 11 χ 13 (t * χ 12 * * * * -(t -t k ) χ 33 (t -t k )y 3	   .

  ).

	It is readily checked that
	lim h→0 +	 sup	ψ(t + h) -ψ(t)

h  ≤ (-2δ + δ 1 )ψ(t).

  1 )t , t ≥ 0. , consider the case where δ 1 = 2δ + 2η. If the LMIs (

	 l	(B.4) e NowThis, together with the following inequality
	0

2ηt [p 1 u 2 (x, t) + p 3 u 2 x (x, t)]dx ≤ V 1 (t) ≤  l 0 [p 1 u 2 (x, 0) + p 3 u 2 x (x, 0)]dx  e (-

2δ+δ

1 )t , t ≥ 0 yields (21).

This work was partially supported by Israel Science Foundation (grant no. 1128/14).


Jonàs Martínez

Random Auxetic Porous Materials from Parametric Growth Processes

Keywords: Porous media, inverse design, random auxetics, growth processes

We introduce a computational approach to optimize random porous materials via parametric growth processes. We focus on the problem of minimizing the Poisson's ratio of a two-phase porous random material, which results in an auxetic material. Initially, we perform a parametric optimization of the growth process. Afterward, the optimized parametric growth process implicitly generates an auxetic random material. Namely, the growth process intrinsically entails the formation of an auxetic material. Our approach enables the computation of large-scale auxetic random materials in commodity computers. We also provide numerical results indicating that the computed auxetic materials have close to isotropic linear elastic behavior and physical tests revealing auxetic behavior.

Introduction

Auxetic structures are materials with a negative Poisson's ratio: when stretched, they expand perpendicularly to the applied force [START_REF] Kolpakov | Determination of the average characteristics of elastic frameworks[END_REF][START_REF] Lakes | Foam structures with a negative Poisson's ratio[END_REF], a seemingly counterintuitive property. Auxetic materials find applications in multiple fields thanks to their excellent shock absorption, fracture toughness, or vibrational absorption [START_REF] Yang | Review on auxetic materials[END_REF][START_REF] Saxena | Three decades of auxetics research -materials with negative Poisson's ratio: a review[END_REF][START_REF] Ma Kolken | Auxetic mechanical metamaterials[END_REF][START_REF] Lakes | Negative-Poisson's-ratio materials: auxetic solids[END_REF][START_REF] Ren | Auxetic metamaterials and structures: A review[END_REF][START_REF] Qin | Negative Poisson's ratio in two-dimensional honeycomb structures[END_REF], among other reasons. A substantial amount of research is devoted to designing auxetic mechanical materials [START_REF] Ma Kolken | Auxetic mechanical metamaterials[END_REF][START_REF] Dagdelen | Computational prediction of new auxetic materials[END_REF][START_REF] Vogiatzis | Topology optimization of multi-material negative poisson's ratio metamaterials using a reconciled level set method[END_REF] that derive their physical properties from the particular arrangement of their small-scale geometry. Recent manufacturing technologies can fabricate complex small-scale structures and therefore manufacture auxetic materials.

Random materials offer some significant advantages compared to the more widespread periodic materials. In particular, they are more resilient to fabrication-related symmetry-breaking imperfections [START_REF] Carlos M Portela | Extreme mechanical resilience of self-assembled nanolabyrinthine materials[END_REF], can smoothly and seamlessly grade material properties [START_REF] Kumar | Inverse-designed spinodoid metamaterials[END_REF], are well suited to manufacture isotropic structures [START_REF] Mullen | Selective laser melting: A unit cell approach for the manufacture of porous, titanium, bone in-growth constructs, suitable for orthopedic applications. ii. randomized structures[END_REF][START_REF] Hossain | Controlling and testing anisotropy in additively manufactured stochastic structures[END_REF], are excellent candidates for energy-absorbing applications [START_REF] Da Chen | Dynamic response and energy absorption of functionally graded porous structures[END_REF][START_REF] Mueller | Energy absorption properties of periodic and stochastic 3D lattice materials[END_REF][START_REF] Ibrahim | Post-yield performance of additive manufactured cellular lattice structures[END_REF], and allow to compute the material geometry efficiently [START_REF] Martínez | Procedural Voronoi foams for additive manufacturing[END_REF]. While a repeating periodic structure defines most auxetic materials, a distinct line of research is otherwise interested in random auxetic materials [START_REF] Mj Mirzaali | Auxeticity and stiffness of random networks: Lessons for the rational design of 3D printed mechanical metamaterials[END_REF] since they offer certain advantages over periodic structures [START_REF] Rayneau-Kirkhope | Density scaling in the mechanics of a disordered mechanical meta-material[END_REF][START_REF] Zerhouni | Numerically-aided 3D printed random isotropic porous materials approaching the Hashin-Shtrikman bounds[END_REF][START_REF] Jagdish | Design of polymeric auxetic matrices for improved mechanical coupling in lead-free piezocomposites[END_REF]. Auxetic polymeric foams [START_REF] Lakes | Foam structures with a negative Poisson's ratio[END_REF][START_REF] Bd Caddock | Microporous materials with negative Poisson's ratios. I. Microstructure and mechanical properties[END_REF] were reported in the 80s and are widely used in industrial applications. The most common process to obtain auxetic foams is compressing a conventional flexible cellular foam to force the cell ribs to buckle, producing a reentrant structure heated to its softening temperature [START_REF] Chan | Fabrication methods for auxetic foams[END_REF][START_REF] Alderson | Manufacture and characterisation of thin flat and curved auxetic foam sheets[END_REF]. The geometry of cellular foams is usually idealized and modeled with Voronoi diagrams [START_REF] Lorna | Cellular solids: structure and properties[END_REF], and some works studying auxetic foams start by modeling a compressed Voronoi diagram [START_REF] Li | Negative Poisson's ratio in 2D Voronoi cellular solids by biaxial compression: a numerical study[END_REF][START_REF] Gao | Numerical analysis of the mechanical properties of 3D random Voronoi structures with negative Poisson's ratio[END_REF]. Alternative methods to produce random auxetic materials consider sheets with random perforations [START_REF] Joseph N Grima | Auxetic perforated mechanical metamaterials with randomly oriented cuts[END_REF] or random fiber networks produced by electrospinning [START_REF] Domaschke | Random auxetics from buckling fibre networks[END_REF].

A recent research direction tackles the minimization of Poisson's ratio of finite two-dimensional random networks consisting of nodes connected by bonds. Reid et al. [START_REF] Daniel R Reid | Auxetic metamaterials from disordered networks[END_REF] proposed to prune network bonds iteratively, further improved in [START_REF] Daniel R Reid | Ideal isotropic auxetic networks from random networks[END_REF] by modifying, in addition, the position of the network nodes and the stiffness of the network bonds. Hagh et al. [START_REF] Varda | Disordered auxetic networks with no reentrant polygons[END_REF] starts from a planar triangulation and iteratively removes bonds while avoiding the creation of reentrant polygons. Liu et al. [START_REF] Liu | Realizing negative Poisson's ratio in spring networks with close-packed lattice geometries[END_REF] proposed an iterative change of the bond stiffness. Pashine et al. [START_REF] Pashine | Directed aging, memory, and nature's greed[END_REF] presented a network aging process that decreases the Poisson's ratio. All these methods assume that the network's elastic behavior is close to isotropic during the iterative optimization process.

We advocate for a different computationally efficient perspective on computing two-dimensional auxetic random materials. We propose to compute a porous material [START_REF] Schroeder | Computer-aided design of porous artifacts[END_REF] with a solid and void phase through a parametric growth process. The growth process is compactly defined by a random point process and two parametric functions controlling the growth law. Our goal is to optimize these two functions' parameters to minimize the Poisson's ratio of the resulting porous material. Importantly, this optimization is only carried out once and allows the immediate generation of random auxetic materials afterward through the process of growth, requiring no further optimization process. This is in opposition to prior methods [START_REF] Daniel R Reid | Auxetic metamaterials from disordered networks[END_REF][START_REF] Daniel R Reid | Ideal isotropic auxetic networks from random networks[END_REF][START_REF] Varda | Disordered auxetic networks with no reentrant polygons[END_REF][START_REF] Liu | Realizing negative Poisson's ratio in spring networks with close-packed lattice geometries[END_REF][START_REF] Pashine | Directed aging, memory, and nature's greed[END_REF] that proceed with iterative optimization, rendering the computation infeasible for large-scale materials. We enable the computation in commodity computers of large-scale random auxetic porous materials with close to isotropic behavior (see Figure 1).

Our research is influenced by [START_REF] Martínez | Starshaped metrics for mechanical metamaterial design[END_REF] that considered a parameterized growth process to obtain periodic porous materials. Notwithstanding, a number of differences distinguish us from [START_REF] Martínez | Starshaped metrics for mechanical metamaterial design[END_REF]. First, we consider the case of random porous materials. Second, we formulate a parametric optimization to deal with the particular case of minimizing the Poisson's ratio. Third, we introduce a novel parameterized set that controls the cease of cell growth and considerably improves the results. Finally, we suggest a hypothesis on the link between certain symmetries of the growth process and the emergence of random porous materials close to isotropic elasticity.

Growth process

We consider a growth process [START_REF] Chiu | Stochastic geometry and its applications[END_REF] in which nuclei are born at the same time. Each nucleus, a point x ∈ R 2 , is the origin of a cell. Cells grow according to a growth law and are forbidden to overlap. At the end of the growth, the union of all cells corresponds to the porous material's void phase, while its complement in R 2 corresponds to the solid phase. We start by describing how the nuclei are placed in space and how the law of growth is defined.

A two-dimensional point process Φ ⊂ R 2 defines the nuclei. We consider two different point processes: a homogeneous Poisson point process with intensity λ > 0, and a hardcore point process in which the points are forbidden to lie closer together than a minimum distance D > 0. In particular, we consider the random sequential adsorption (RSA) model [START_REF] Widom | Random sequential addition of hard spheres to a volume[END_REF], where disks with radius D 2 are iteratively added as long as they do not overlap any previously added disk. The center of the disks corresponds to the point process, and its intensity is λ = 4A πD 2 , where A ≈ 0.547 is the area fraction where saturation occurs [START_REF] Chiu | Stochastic geometry and its applications[END_REF].

Cells grow through a uniform scaling of a compact set S ⊂ R 2 centered on each nucleus. The growth ceases whenever and wherever a cell comes into contact with another (see Figure 2). This type of growth process has been well studied in the literature when S is a convex set and is also known as the Voronoi growth model [START_REF] Okabe | Spatial tessellations: concepts and applications of Voronoi diagrams[END_REF]. Our central idea is to allow S to be non-convex, enabling the growth of porous materials with more diverse shapes.

In addition, we control how close to each other cells can grow, with what we call the maximum Euclidean radius of growth (see Figure 3). Let B r ⊂ R 2 be a closed Euclidean disk with radius r and centered at the origin. Let x be a point in the boundary of a growing cell. The growth ceases at x if x + B r intersects another growing cell.

We start by defining S. Let S ⊂ R 2 be a compact star-shaped set with respect to the origin O = (0, 0). For any point x ∈ S, the line segment [x, O] is contained in S. The boundary of S can be parameterized in polar coordinates by a continuous and periodic function ψ S : [0, 2π] → [r min , r max ], where 0 < r min ≤ r max are the minimum and the maximum Euclidean distance from O to the boundary of S.

Following [START_REF] Martínez | Starshaped metrics for mechanical metamaterial design[END_REF], we parameterize S with a few known values (α i , l i ) i>0 of the function ψ S (α i ) = l i , and we call them radial spans. For simplicity, we always consider m > 0 equally spaced radial spans separated by an angle of β m = 2π m and starting at an angle of zero α 1 = 0. Let α jk ∈ [0, 2π] be an angle we seek to interpolate, lying between (α j , l j ) and (α k , l k ), such that j = α jk β m and k = (j + 1) mod m. The relative position t of α jk between α j (t = 0) and α k (t = 1) is given by the fractional part t = α jk β mα jk β m ∈ [0, 1). Then, we interpolate the value of ψ S (α jk ) with the following cubic Hermite spline on the unit interval

ψ S (α jk ) = (2t 3 -3t 2 + 1)l j + (-2t 3 + 3t 2 )l k (1)
The interpolation ensures that ψ S (α jk ) lies between l j and l k by setting the starting and ending tangents of the spline to zero. Besides, we can impose rotational symmetries to the interpolation. Figure 4 illustrates the interpolation procedure. The distance induced by S from a point p to a point x [START_REF] Martínez | Starshaped metrics for mechanical metamaterial design[END_REF] allows us to define the cell growth process compactly

d S (p, x) = x -p ψ S ( (x -p)) (2)
where (•) is the angle of a vector with respect to the horizontal axis, and • is the Euclidean norm. d S (p, x) is the factor by which we have to uniformly scale S, centered on p, such that x lies on its boundary. We always consider that p is a cell nucleus and x is any point of the cell emanating from p. A star-shaped set S * parameterizes the maximum Euclidean radius of growth. We consider a polar function ψ S * : [0, 2π] → [r * min , r * max ], where 0 < r * min ≤ r * max are the minimum and maximum Euclidean radii of growth. Let x be a point in the boundary of a growing cell emanating from p ∈ Φ. The growth ceases at x if x + B ψ S * ( (x-p)) intersects any other growing cell. Thus, S * determines when to cease one cell's growth process as it approaches other cells (see Figure 3).

The minimal distance between two points in a Poisson point process can be arbitrarily small. Thus, a cell can be empty as its growth was forbidden through S * from the beginning of the growth process. For a hardcore point process, the minimal distance between two points is at least D. Thus, it suffices to enforce D > r * max t in order to guarantee that all cells are allowed to grow.

It is challenging to formulate a continuous growth law [3] since both S and S * are not necessarily convex, and S * imposes an atypical constraint on the growth process. We instead resort to a discrete growth process formulation [START_REF] Aletti | Statistical aspects of fuzzy monotone set-valued stochastic processes. application to birth-and-growth processes[END_REF], where the cells grow in discrete steps. Therefore, our growth process produces a discrete random set. The general principle is to assign each cell a disjoint discrete set of points. More precisely, we consider an integer lattice

L W = Z 2 ∩ W where W is a bounded subset of R 2 . Let Φ W = Φ∩W = (p 1 , .
. . , p g ) be the subset of nuclei in W . The set of lattice points belonging to the cells is given by the union

C W = g i=1 C i W ⊂ L W ,
where C i W corresponds to the lattice cell points associated with the nuclei p i . Lattice points are incrementally added to C W throughout the discrete growth process. Let S 1 ⊂ R 2 be the unit square centered at the origin and let ⊕ denote the Minkowski sum operator.

We consider that the porous material's void phase is given by C W ⊕ S 1 , while the solid one is given by (L W \ C W ) ⊕ S 1 . If W is a rectangle, both phases can be compactly represented with a binary digital image. In all illus-trations, the void phase is colored in white and the solid phase in black. Figures 5 and6 show some examples of the growth process. The size of the solid phase can be regulated by S * . For instance, lower values of r * min will lead to narrower porous materials.

S

Cell growth Porous material In the following, we describe in more detail Algorithm 1, which defines the process of discrete growth. An ordered set Q contains tuples {x, p} where the first element is a lattice point of x ∈ L W and the second one is a nucleus point of p ∈ Φ W . Q is ordered according to the distance d S (p, x). The tuples in Q identify the points where cell growth might take place. The growth process is divided into two parts, the nucleation phase and the growth phase.

Nucleation. Before the growth starts, for all p i ∈ Φ W , we insert in Q the tuple {x i , p i }, where x i ∈ L W is the closest lattice point to p i . Ideally, L W should be sufficiently large, so no two points in Φ W are closest to the same lattice point in L W .

Growth. The discrete growth process is simulated as follows. Let {x, p i } be the tuple in Q with the smallest distance d S (p i , x) among all tuples.

-(Line 7) First, we remove {x, p i } from Q.

-(Lines 8-11) Second, if x does not already belong to any cell and the Euclidean disk B ψ S * ( (x-pi)) does not contain any other cell point, we add x in C i W . That is, the lattice point x now belongs to the cell emanating from p i .

-(Lines 12-14) Finally, candidate points of growth around a local neighborhood of p are inserted in Q.

More precisely, we insert in Q the four tuples {x n , p i }, where x n is a lattice point in the 4-connected neighborhood of x in L W .

The previous steps are iteratively done until Q is empty.

Algorithm 1 discreteGrowthProcess(L W , Φ W , S, S * ) 1: C i W ← ∅ for all C i W ∈ C W 2: for p i ∈ Φ W do Nucleation 3:
Let x i ∈ L W be the closest lattice point to p i 4:

Insert in Q the tuple {x i , p i } 5: end for

6: while Q = ∅ do Growth 7: Let {x, p i } ∈ Q with smallest distance d S (p i , x) 8:
Remove {x,

p i } from Q 9: if x / ∈ C W then 10: N ← L W ∩ x + B ψ S * ( (x-pi))
11:

if (y / ∈ C W ) ∨ (y ∈ C i W )
for all y ∈ N then 12:

C i W ← C i W ∪ {x} 13:
for x n ∈ neighborhood of x in L W do end if 20: end while Computational complexity of Algorithm 1. Let n be the number of points of L W and g be the number of nuclei of Φ W . Since we assume that n > g, the number of iterations of the nucleation loop (lines 2-5) and the growth loop (lines 6-20) is of order n. The size of Q is O(n). If Q is implemented as a heap priority queue, then the insertion and removal of elements in Q can be done in logarithmic time (lines 4, 7, and 15). The loop of line 10 (i.e., for all y ∈ N ) has O(1) time complexity since it depends on the constant size of S * . The loop of line 13 also has O(1) time complexity (four iterations). Thus, Algorithm 1 has O(n log(n)) worst-case time complexity. In practice, we observe a quasi-linear runtime with respect to n (see Section 4.3) primarily because the size of the queue Q remains comparatively small in proportion to n. Periodic boundary conditions. Optionally, when W is a rectangle, periodic boundary conditions can be easily imposed by considering that L W is a periodic rectangular lattice as well. If Φ is an RSA point process, the random adsorption of disks must be performed considering the periodicity. All results shown in this article feature periodic boundary conditions. Cell regularization. Algorithm 1 is guaranteed to produce cells C i W that are 4-connected sets. However, we have observed that when S * is not convex, a cell may have a genus higher than zero (i.e., has "holes"). Holes lead to a solid phase composed of more than one connected component, which is undesirable. To resolve this, we fill any hole that a cell may have. Formally, we consider the complement of the unbounded component of the complement of a cell (see Figure 7). We denominate this last step cell regularization. Numerical results show that cell regularization has, on average, little impact (see Section 4.1).

S S *

Cell growth Regularized

Elastic behavior of random porous materials

We study the elastic behavior, assuming linear elasticity and plane stress [START_REF] Thorpe | New results in the theory of elasticity for two-dimensional composites[END_REF], considering that the stresses in the orthogonal direction to the plane are negligibly small (e.g., a thin plate). We apply a procedure known as stochastic homogenization [START_REF] Blanc | Some variance reduction methods for numerical stochastic homogenization[END_REF] that seeks to approximate the average elastic behavior of random materials. Hooke's law represents the elastic behavior of materials with the linear relation σ = C , where σ is the stress, is the strain, and C is the elasticity tensor (a symmetric and positive definite matrix), in Voigt notation [START_REF] Voigt | Lehrbuch der kristallphysik[END_REF]

  σ 1 σ 2 σ 12   =   c 11 c 12 c 13 c 12 c 22 c 23 c 13 c 23 c 33     1 2 2 12   (3)
Periodic homogenization [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] seeks to find the elasticity tensor characterizing a periodic composite material defined from a periodic cell, as the length of the cell tends to zero. In the periodic case, the homogenization is defined on a finite domain (periodic cell), while in the random case is defined on the whole space R 2 and cannot be reduced to a problem posed on a finite domain [START_REF] Bourgeat | Approximations of effective coefficients in stochastic homogenization[END_REF].

We approximate the homogenized elasticity tensor coefficients of a random porous material using the so-called "cut-off" techniques. Consider a periodized cut [0, s] 2 (a square) of a stationary random material and its corresponding homogenized elasticity tensor C s . It was shown in [START_REF] Bourgeat | Approximations of effective coefficients in stochastic homogenization[END_REF] that lim s→∞ C s converges, almost surely. Nevertheless, s needs to be bounded for computation. In our discrete setting, the cut size s is a positive integer that corresponds to the pixel resolution s × s of the porous material image. In all results, we consider a constant D = 80 for the RSA point process and λ = 4A πD 2 ≈ 0.000 11 for the Poisson point process. Therefore, both point processes have similar intensities.

We compute the homogenized elasticity tensor C s with a publicly available numerical homogenization method based on the finite element method and considering plane stress properties [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF]. The input is an image with values identifying each distinct phase of the periodic material. We impose periodic boundary conditions on the growth process to produce a suitable input image. The solid phase is a linear elastic isotropic material with unit Young's modulus E = 1.0 and Poisson's ratio v = 0.3.

Let E[•] be the expectation of a random variable. We estimate the value of E[C s ] with the average of k independent realizations of the random porous material giving

(C s ) i , 1 < i ≤ k [6] E[C s ] ≈ C s = 1 k k i=1 (C s ) i (4)
Let std(C s ) be the standard deviation of the k independent realizations. A range of plausible values for E[C s ] is

C s ± 1.96 std(C s ) √ k
for a confidence interval with probability 95% (see Figure 8). In order to evaluate the accuracy of the stochastic homogenization [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF], we examine the statistical fluctuations of the material porosity p ∈ [0, 1] (the area fraction of the void phase) for a given k and s. In particular, we consider the coefficient of variation c v (p) = std(p) p that indicates the extent of variability of the porosity in relation to the mean. Increasing the value of k and s most likely decreases c v (p). Following [START_REF] Stroeven | Numerical determination of representative volumes for granular materials[END_REF], we consider that the approximation is sufficiently accurate if c v (p) is below a maximum predefined one. In all the following results, we always consider k = 60 realizations.

Isotropic elasticity

Since our goal is to minimize the Poisson's ratio of an isotropic material, we give a brief presentation of how we determine if the average tensor C s is close to isotropic, and if so, what are its defining parameters E and v. An isotropic material has an elasticity tensor C iso with only two independent components, c iso 11 and c iso

For plane elasticity, the relation between the tensor components c iso 11 , c iso 12 and the Young's moduli E > 0 and the Poisson's ratio v ∈ [-1, 1] is given by

-→ i = 1 i = 2 i = 3 i = 4 i =
E = (c iso 11 ) 2 -(c iso 12 ) 2 c iso 11 v = c iso 12 c iso 11 (6)
In practice, it never occurs that the average tensor C s is perfectly isotropic. Thus, we seek to determine its closest isotropic tensor under some norm [START_REF] Moakher | The closest elastic tensor of arbitrary symmetry to an elasticity tensor of lower symmetry[END_REF]. We consider the Frobenius norm for a 3 × 3 matrix A

A F = 3 i=1 3 j=1 |a ij | 2 (7)
and minimize the Frobenius norm between C iso and the average tensor

C s min c iso 11 ,c iso 12 C s -C iso F (8)
We provide in the Appendix A the closed-form expression of c iso 11 , c iso 12 from C iso minimizing Equation 8. We also need a measure of the deviation of C s from perfect isotropy. We consider the normalized measure δ iso used in practice [START_REF] Vp Tran | Stochastic modeling of mesoscopic elasticity random field[END_REF][START_REF] Zerhouni | Numerically-aided 3D printed random isotropic porous materials approaching the Hashin-Shtrikman bounds[END_REF]

δ iso = C s -C iso F C s F ≥ 0 ( 9
)
where δ iso = 0 indicates perfect isotropy and increasing values indicate divergence from isotropy.

Results

Three-fold symmetry and isotropy

It is known that two-dimensional periodic structures with three-fold rotational symmetry (invariant to rotations of 2π

3 ) lead to an isotropic linear elastic behavior [START_REF] Forte | Symmetry classes for elasticity tensors[END_REF]. We hypothesize that a resembling property holds in our instance of random growth.

A two-dimensional point process Φ is stationary if Φ+x has the same probability distribution as Φ for all x ∈ R 2 . Φ is isotropic if R(a)Φ has the same probability distribution as Φ, for all rotations R(a) of angle a ∈ R around the origin. The Poisson point process and the RSA point process are stationary and isotropic point processes [START_REF] Chiu | Stochastic geometry and its applications[END_REF][START_REF] Stoyan | Random sequential adsorption: relationship to dead leaves and characterization of variability[END_REF].

Let us consider a rotation of R( 2π 3 n), for n ∈ Z. The probability distribution of Φ is invariant to a rotation R( 2π 3 n). If S and S * are three-fold symmetric, then R( 2π 3 n)S = S (same for S * ). Given a point p ∈ Φ then we have

R( 2π 3 n)(p + S) = R( 2π 3 n)p + R( 2π 3 n)S = R( 2π 3 n)p + S (10
) i.e., the rotation of S centered at p around the origin keeps S and S * invariant. This suggests that the probability distribution of the discrete random set given by the growth process is close to being stationary and invariant to rotations R( 2π3 n), which leads us to the following hypothesis. Hypothesis 1. Let Φ be a stationary and isotropic point process. Let S and S * be three-fold rotational symmetric. For s → ∞, the discrete growth process described in this article gives a porous material with δ iso close to zero. Namely, the porous material has close to isotropic linear elastic behavior.

We numerically tested Hypothesis 1 for an extensive set of porous materials with different sets of S and S * , and increasing size s (see Figure 9). There is strong evidence that as s increases, the average value of δ iso decreases, and it appears to tend towards a value close to zero.

In all the subsequent results, we consider a square image of side length s = 1200 pixels. This size gave among all the tests performed in Figure 9 a maximum coefficient of variation of porosity of c v (p) = 0.035 19 for the Poisson point process, and of c v (p) = 0.011 18 for the RSA point process.

We also evaluate the impact of cell regularization over the porosity by measuring the relative number of filled cell holes. Let p be the porosity of the porous material before performing the cell regularization. We consider the relative measure ζ p, p = p-p p ≥ 0. For all the results in Figure 9 the maximum value of ζ p, p is 0.061 57 and the mean is 0.004 33 for the Poisson point process, and maximum 0.019 48 and mean 0.000 41 for the RSA point process. This indicates that the effect of cell regularization on the porosity is almost negligible, particularly for the RSA point process.

Parametric optimization

By this point, we have all the ingredients to formulate the parametric optimization of the growth process. Let C s (ψ S , ψ S * ) be the average elasticity tensor of a porous material parameterized by the functions ψ S (star-shaped distance) and ψ S * (maximum Euclidean radius of growth). Let F : M → R be the objective function, where M is a 3 × 3 symmetric matrix. Our objective is to minimize arg min

ψ S ,ψ S * F(C s (ψ S , ψ S * )) = c iso 12 c iso 11 (11)
where c iso 11 , c iso 12 are derived from C s (ψ S , ψ S * ) considering Equation A.4. To attain a physical behavior close to isotropic elasticity, we constrain S and S * to be three-fold rotational symmetric (see Hypothesis 1).

The input of function F is expensive to evaluate and has no closed-form. In this case, Bayesian optimization [START_REF] Močkus | On Bayesian methods for seeking the extremum[END_REF] is often a good choice to tackle the optimization. In particular, we have used Bayesian optimization using Gaussian processes as implemented in the Python library scikitoptimize [START_REF] Head Tim | scikitoptimize 0[END_REF]. Both functions ψ S and ψ S * are parameterized with m = 8 radial spans. We set the first radial span of ψ S to a constant value of 1, removing an unnecessary degree of freedom. Finally, both ψ S and ψ S * are concatenated into a single optimization variable to optimize ψ S and ψ S * simultaneously. We set a maximum number of 150 evaluations of function F and 50 random initial guesses.

Figures 10 to 12 provide minimization results of the Poisson's ratio with increasingly larger optimization domains of ψ S * and a constant optimization domain of ψ S ∈ [0.05, 1.0]. We distinguish two clear trends from the optimization results. First, the RSA point process can achieve a lower Poisson's ratio v compared to the Poisson point process. Second, increasing the function ψ S * domain leads to a lower Poisson's ratio since we expand the space of achievable porous materials. Overall, we recommend using a hardcore point process instead of a Poisson point process to reach a lower Poisson's ratio and lower variance of the material properties using the same area (see Figure 9).

Computational efficiency

We evaluate the efficiency of the growth process computation by performing some timing tests shown in Table 1. For a range of n = [1, 100] million pixels, we observe a nearly linear correlation between n and computing time. Based on the statistics shown in Table 1, we consider that this is mainly explained because the relative average size of the queue Q is below 4% (column c ○), amortizing the logarithmic cost of insertion and deletion of elements in Q. In addition, the relative time spent computing the cell regularization is below 1% (column a ○), showing that the bulk of the computational cost is in the growth algorithm.

Finally, we observe different scaling laws among the tests, depending on the input sets S and S * and the point process type. We infer that these differences are explained mainly by a trade-off between the number of iterations in the main loop (column b ○) and the average number of iterations in line 11 of Algorithm 1 (column d ○). Both statistics are tightly related to S * : a larger area of S * implies fewer iterations of b ○ at the expense of more iterations of d ○.

RSA point process

Poisson point process

RSA Poisson ψ S * a ○ b ○ c ○ d ○ a ○ b ○ c ○ d ○ 4.0 0.995% 453%

Conclusions and future work

We have presented a computational approach based on parametric growth processes to minimize the Poisson's ratio of random porous materials. The results attain a negative Poisson's ratio while providing a numerical analysis showing close to linear elastic isotropic behavior. Increasing the number of parameters of ψ S and ψ S * and their domain further improves the minimization. In practice, our approach can compute large-scale porous materials in a reasonable time. This capability is most valuable for applications that need to generate large quantities of random porous materials.

It is possible to generate heterogeneous materials with a parametric growth process by considering that both functions ψ S and ψ S * vary in the space. One of our approach's significant advantages over periodic structures is that the transitions between different material properties do not need to be handled [START_REF] Garner | Compatibility in microstructural optimization for additive manufacturing[END_REF][START_REF] Hu | Cellular structure design based on free material optimization under connectivity control[END_REF], as they are intrinsically cap-tured by the parameters governing the growth process. The gradations can be implicitly done either by constant (see Figure 15) or linear interpolation (see Figure 16) of the star-shaped sets. Figure 17 shows a 3D rendering of thin sheets obtained by our contouring method (see Section 4.4).

Manufacturing processes impose morphological constraints. Having a minimum local size of the solid phase is a common one. For a varying function ψ S * , we may obtain porous materials whose solid phase has features with minimal local size (i.e., dimension of one-pixel, see Figure 14). Thus, it would be helpful to study how to ensure a minimum local size of the solid phase everywhere, either as a postprocessing step or as part of the growth process.

Our approach also opens up different opportunities for future work. First, to parallelize Algorithm 1 and significantly increase its performance. Second, to explore the much wider space of porous materials given by S and S * having other symmetries and other point distributions than those considered in this article. Third, to extend Algorithm 1 to the 3D case by considering a 3D growth process and adapting the parameterization of S and S * . The inverse design of random 3D porous materials (e.g., random foams) is of particular interest for acoustic absorption [START_REF] Han | Acoustic absorption behaviour of an open-celled aluminium foam[END_REF] or tissue engineering [START_REF] Lanza | Principles of tissue engineering[END_REF]. Finally, it would be interesting to perform further physical tests of the results, considering larger deformations of the porous material (non-linear regime) and mechanical properties such as stress failure or buckling.

Our approach may also open the door to other applications besides minimizing the Poisson's ratio. For example, to tackle the inverse design and optimization of porous materials [START_REF] Ronellenfitsch | Inverse design of discrete mechanical metamaterials[END_REF][START_REF] Kumar | Inverse-designed spinodoid metamaterials[END_REF] using our growth process. We conjecture that it will be sufficient to replace the objective function of Equation 11 with another one of interest in most cases.

Code and data availability

An open-source implementation of our approach can be found at1 . We provide a script to automate the generation of all paper results and figures, apart from Figure 3 for each sheet. The red dotted line is superimposed on the image to appreciate the lateral expansion better. Since our numerical results consider linear elasticity, we apply a relatively small uniform stretching of 10 mm (stretch ratio 5%) to all the sheets. The Poisson's ratio v corresponds to the mean numerical value obtained with stochastic homogenization. We visually observe a monotonic correlation between the lateral displacement and the corresponding numerical value of v.

  v

Figure 1 :Figure 2 :

 12 Figure 1: From left to right: large-scale random auxetic porous material with increasing close-up views. The growth process has 43590 nuclei, and the computation took 26 minutes on a laptop using a single CPU core.

Figure 3 :Figure 4 :

 34 Figure 3: Growth of two cells. The growth ceases at x since the disk x + Br intersects another growing cell.

Figure 5 :Figure 6 :

 56 Figure5: The growth process is given by the point process of Figure2and the uniform scaling of S at the left. In this example, S * is a Euclidean disk.

Figure 7 :

 7 Figure 7: Left: S and S * governing the growth process of two nuclei. Middle: evolution of the growth. In this example, the holes arise because of the high anisotropy of S and S * . Right: cell regularization, the two cell holes are filled.

33 Figure 8 :

 338 Figure8: Approximation of the elasticity tensor of the random porous material of Figure6in a square with size s = 1200, and 120 independent realizations. Top: Some realizations. Bottom: mean and confidence interval of each tensor component. As the number of realizations i increases, the mean stabilizes, and the confidence interval becomes narrow. The coefficient of variation of porosity is cv(p) = 0.004 09.

Figure 9 :

 9 Figure 9: Analysis of the deviation of isotropy δ iso and coefficient of variation of the porosity cv(p) for a large number of porous materials and two different point processes (left and right columns). S and S * are three-fold rotational symmetric. The colored boxes encompass the interquartile range of values of δ iso . We consider 4 different sizes s, going from 300 to 1200 For each point process and size, we enumerate 300 different random radial spans in the range ψ S ∈ [0.05, 1.0] and ψ S * ∈ [4.0, 40.0], and 8 different equally spaced angles in [0, 2π 3 ] (three-fold symmetry) considering a total of 2 • 4 • 300 = 2400 different porous materials. We observe a clear trend of the average value of δ iso and cv(p) getting closer to zero as the size s of the domain increases. We also notice lower average values with the RSA point process compared to the Poisson point process.

Figure 10 :Figure 11 :

 1011 Figure 10: Optimization with constant ψ S * = 4.0.

Figure 12 :

 12 Figure 12: Optimization with variable ψ S * ∈ [4.0, 40.0]

  , we perform 10 tests with increasing size n and plot the corresponding time in minutes. The table below shows the following additional statistics (considering the average value between all the tests). a ○: Relative time spent computing the regularization with respect to the total time of computation. b ○: Number of iterations in the main loop (line 6 of Algorithm 1) divided by n. c ○: Average size of the queue Q divided by n. d ○: Average number of iterations in line 11 of Algorithm 1.

ψFigure 13 :

 13 Figure13: Visualization of the Poisson's ratio of four laser-cut rubber sheets. We show the picture in rest (left) and stretched position (right) for each sheet. The red dotted line is superimposed on the image to appreciate the lateral expansion better. Since our numerical results consider linear elasticity, we apply a relatively small uniform stretching of 10 mm (stretch ratio 5%) to all the sheets. The Poisson's ratio v corresponds to the mean numerical value obtained with stochastic homogenization. We visually observe a monotonic correlation between the lateral displacement and the corresponding numerical value of v.

Figure 14 :

 14 Figure 14: The solid phase in black color has a little feature inside the red square (left) that is visually magnified (right).

Figure 15 :

 15 Figure 15: Spatially varying growth process with two different growth functions. Left: Input heterogeneous field. The darker regions correspond to the optimized ψ S in Figure 10, while the rest correspond to growth according to the Euclidean disk.

Figure 16 :

 16 Figure 16: Spatially varying growth process with three different growth functions linearly interpolated. Left: Input field α ∈ [0, 2] to control the linear interpolation between the optimized sets from Figure 10 (S α=2 ), Figure 11 (S α=1 ), and Figure 12 (S α=0 ). More precisely, for α ∈ [0, 1) we consider the interpolated set α(S α=1 ) + (1 -α)(S α=2 ), and for α ∈ [1, 2] we consider (α -1)(S α=0 ) + (2 -α)(S α=1 ).

Figure 17 :

 17 Figure 17: 3D rendering as thin sheet plates of the results of Figure 15 (left) and Figure 16 (right).

Table 1 :

 1 Timing tests were performed on images with sizes ranging n = [1, 100] million pixels. We compute the growth process with a laptop equipped with an Intel i7-8650U. For each optimization result in Figures 10 to 12

	3.07% 79.4 0.998% 456% 3.19% 79.2
	[4.0, 15.0] 0.367% 342% 2.98% 314 0.402% 355% 2.59% 294
	[4.0, 40.0] 0.363% 226% 2.52% 759 0.171% 226% 1.47% 2930

https://github.com/mfx-inria/auxeticgrowthprocess2d

Acknowledgments

This work was partly supported by ANR MuFFin (ANR-17-CE10-0002). Thanks to Pierre Bedell and Pierre-Alexandre Hugron for their help in doing the physical experimental tests.

Appendix A. Derivation of c iso 11 and c iso

For the three-dimensional case, the minimizing values of c iso 11 , c iso 12 were derived, for instance, in [START_REF] Norris | The isotropic material closest to a given anisotropic material[END_REF]. We were not able to find any previous work concerning the more straightforward two-dimensional case. Nevertheless, the values minimizing Equation 8 can be derived as follows. Consider the function f with two arguments c iso 11 and c iso 12 (A.1) The Hessian matrix

is positive definite everywhere, and it follows that f is a convex function. Since found by solving the system of linear equations from Equation A.2 being equal to zero.