Pierre Béjot
email: pierre.bejot@u-bourgogne.fr

Bertrand Kibler

Quadrics for structuring invariant space-time wavepackets

Keywords: Space-time wavepackets, orbital angular momentum, optical waveguides, light structuring

Space-time light structuring has emerged as a very powerful tool for controlling the propagation dynamics of pulsed beam. The ability to manipulate and generate space-time distributions of light has been remarkably enhanced in past few years, letting envision applications across the entire spectrum of optics. Space-time invariant optical wavepackets manipulated up to now are usually twodimensional objects (one space dimension and time) whose mode-resolved spectra lie in a conical section. Using simple symmetry and invariance principles, we show that such wavepackets are particular cases of more general three-dimensional structures whose space-time frequencies lie on quadric surfaces. Our proposed framework allows here classifying invariant space-time wavepackets localized in all dimensions, in any group-velocity dispersion regime, both in bulk and waveguides. Particular emphasis is placed on longitudinal orbital angular momentum-carrying space-time wavepackets. This unprecedented theoretical approach opens the way for versatile synthesizing of space-time optics.

Introduction. Symmetry and invariant principles play a fundamental role in modern physics. [START_REF] Wigner | Symmetry and conservation laws[END_REF][START_REF] Gross | The role of symmetry in fundamental physics[END_REF] Major examples from classical mechanics are the total energy, and linear and angular momenta obeying respective conservation law associated with continuous symmetries (time translation, space translation and rotation). The group theory is the prominent mathematical approach for studying symmetries in various branches ranging from solid state physics to quantum physics. It allows to analyze symmetries and their implications in a unified way for distinct physical systems, including solving or simplifying various problems, such as obtaining analytic solutions of differential equations (e.g. Lie's method). Application of symmetries in optical physics can be easily found in geometrical optics and optical aberrations, [START_REF] Khan | Group theory in geometrical optics[END_REF][START_REF] Dragt | Lie algebraic theory of geometrical optics and optical aberrations[END_REF] as well as the derivation of non-diffracting solutions in free space. [START_REF] Boyer | Symmetry and separation of variables for the Helmholtz and Laplace equations[END_REF][START_REF] Wünsche | Generalized Gaussian beam solutions of paraxial optics and their connection to a hidden symmetry[END_REF][START_REF] Rodríguez-Lara | Symmetry in optics and photonics: a group theory approach[END_REF] For the latter, the symmetry group associated with Helmholtz equation and the coordinate systems in which variables separate give rise to various solution families of invariant or accelerating waves and beams such as plane waves, Bessel and Mathieu beams, and Weber waves. [START_REF] Durnin | Exact solutions for nondiffracting beams. I. The scalar theory[END_REF][START_REF] Gutiérrez-Vega | Alternative formulation for invariant optical fields: Mathieu beams[END_REF][START_REF] Bandres | Nondiffracting accelerating waves: Weber waves and parabolic momentum[END_REF] Orbital angular momentum (OAM)carrying beams (i.e., vortex beams) are other possible monochromatic solutions (not always diffraction-free) associated to screw axis symmetry. [START_REF] Shen | Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities[END_REF] Besides the intense research activities on optical vortex beams carrying longitudinal OAM (i.e., the spiral phase twists in the spatial domain), the recent studies on spatiotemporal optical vortices (STOVs) carrying transverse or tilted OAM (here the phase rotates in the spatiotemporal domain) have opened up new opportunities for structuring light. [START_REF] Bliokh | Spatiotemporal vortex beams and angular momentum[END_REF][START_REF] Chong | Generation of spatiotemporal optical vortices with controllable transverse orbital angular momentum[END_REF][START_REF] Hancock | Free-space propagation of spatiotemporal optical vortices[END_REF][START_REF] Hancock | Mode structure and orbital angular momentum of spatiotemporal optical vortex pulses[END_REF][START_REF] Cao | Non-spreading Bessel spatiotemporal optical vortices[END_REF][START_REF] Wan | Photonic orbital angular momentum with controllable orientation[END_REF] More generally, structured light, particle beams, and invariant wavepackets are nowadays indispensable tools in numerous fields of science and technology, [START_REF] Rubinsztein-Dunlop | Roadmap on structured light[END_REF][START_REF] Forbes | Structured light[END_REF][START_REF] Piccardo | Roadmap on multimode light shaping[END_REF][START_REF] Lloyd | Electron vortices: Beams with orbital angular momentum[END_REF][START_REF] Urrutia | Helicon waves in uniform plasmas. IV. Bessel beams, Gendrin beams, and helicons[END_REF][START_REF] Bliokh | Spin and orbital angular momenta of acoustic beams[END_REF][START_REF] Vanacore | Ultrafast generation and control of an electron vortex beam via chiral plasmonic near fields[END_REF] and a new step forward in the advanced control of wave propagation would facilitate the emergence of new functionalities. While most of the above examples have been deeply studied in the monochromatic regime, strong interest is currently shown for space-time (ST) invariant wavepackets. Among them, one can cite space-time light sheets, [START_REF] Kondakci | Diffractionfree space-time light sheets[END_REF] needle pulses, 26 X-and O-waves. [START_REF] Kibler | Discretized conical waves in multimode optical fibers[END_REF][START_REF] Faccio | Conical Waves, Filaments and Nonlinear Filamentation Optics[END_REF][START_REF] Hernandez-Figueroa | Localized Waves[END_REF] Invariant STOVs have also been observed in the anomalous dispersion regime. [START_REF] Cao | Non-spreading Bessel spatiotemporal optical vortices[END_REF] Structuring the space-time characteristics of a field turns out to be a very versatile way to engineer the pulse propagation dynamics, allowing, for instance, to control the group velocity [START_REF] Kondakci | Diffractionfree space-time light sheets[END_REF][START_REF] Kondakci | Optical spacetime wave packets having arbitrary group velocities in free space[END_REF] or group velocity dispersion, [START_REF] Yessenov | Engineering the Optical Vacuum: Arbitrary Magnitude, Sign, and Order of Dispersion in Free Space Using Space-Time Wave Packets[END_REF][START_REF] Hall | Canceling and inverting normal and anomalous group-velocity dispersion using space-time wave packets[END_REF] as well as OAM features [START_REF] Cao | Non-spreading Bessel spatiotemporal optical vortices[END_REF][START_REF] Yessenov | Spacetime wave packets localized in all dimensions[END_REF] of the structured wavepackets. Not only used for sculpting "at will" the propagation dynamics, the framework of invariant ST wavepackets has provided a deeper understanding of the mechanisms responsible for supercontinua generation during the nonlinear propagation of ultrashort pulses. [START_REF] Kibler | Discretized conical waves in multimode optical fibers[END_REF][START_REF] Kolesik | Interpretation of the spectrally resolved far field of femtosecond pulses propagating in bulk nonlinear dispersive media[END_REF][START_REF] Béjot | Spatiotemporal Helicon Wavepackets[END_REF] In a geometrical point of view, all the aforementioned studied space-time structures are two-dimensional and share the particularity to own a frequency-resolved modal decomposition that lies on a conic section. Since extremely fast progresses are currently made for generating more and more complex three-dimensional light structures, [START_REF] Yessenov | Spacetime wave packets localized in all dimensions[END_REF][START_REF] Guo | Structured 3D linear space-time light bullets by nonlocal nanophotonics[END_REF] it then seems particularly important at this point to provide an exhaustive classification of three-dimensional invariant ST wavepackets.

In the present work, we provide a general approach for structuring three-dimensional light ultrashort pulses in both bulk and structured dispersive media. We show that, in two different space coordinates systems (Cartesian and cylindrical), ST wavepackets that are invariant under uniform rectilinear-motion and spatial rotation operations can be built from the coherent superposition of modes that lie on quadric surfaces in the frequency-resolved modal space. The nature of surfaces depends on the considered dispersion regime and the group-velocity of the wavepacket. The geometrical study of the considered quadrics then allows constructing invariant ST wavepackets propagating with arbitrarily chosen group velocity and rotation. Every invariant ST wavepackets studied so far could be investigated within a similar framework. Note that the same analysis could be performed in spherical coordinates [START_REF] Mills | Localized waves with spherical harmonic symmetries[END_REF] or space-time cylindrical coordinates, [START_REF] Cao | Non-spreading Bessel spatiotemporal optical vortices[END_REF] the quest for invariant wavepackets being intertwined with symmetries of the considered medium of propagation.

Propagation equation in the modal representation.

The partial differential equation driving the linear propagation in bulk media of an electric field is given by: ,

where , is the light velocity in vacuum, is the refractive index at the angular frequency , is the transverse Laplacian operator, is the transverse position vector, and is the Fourier transform of with respect to the time coordinate. The very first step for studying z-invariant ST wavepackets is to express the eigenvectors (called hereafter "modes") of the operator . In the Cartesian (resp. cylindrical) case, they simply correspond to monochromatic plane waves [resp. orbital angular momentum Bessel monochromatic waves ] with eigenvalues , where . The evolution of a field along the direction is given by: ,

where is the decomposition coefficients of the field in the basis composed of the eigenvector of Any electric field can be then represented in the modal basis in Cartesian coordinates, and in cylindrical coordinates.

Galilean-transformed (uniform rectilinear motion) invariant wavepackets.

The first example of invariant ST wavepackets are those propagating without any deformation in a frame propagating at their group velocity . In this case, the transformation linking the field at a given propagation distance with the initial field is simply . Finding the spatiotemporal wavepackets invariant under uniform rectilinear motion (and invariant under phase transformation), amounts to find the fields that belong to the kernel of the differential operator , where 1/ is the group velocity of the built wavepacket and is its propagation constant. It immediately follows that the decomposition of rectilinear motion-invariant wavepackets only embeds eigenvectors whose eigenvalues respect: .

(

) 3
One can now study how the mode-resolved spectrum of such invariant fields are structured in bulk media. For simplicity, we will consider the dispersive properties of the medium by means of the usual Taylor expansion around a given arbitrarily chosen frequency and up to the second-order dispersion coefficient. Now, Eq. 3 reads: ,

The above algebraic expression in 3D-space refers to a quadric surface of dimension 2 (note that higher order terms in have been neglected). Such quadrics are easily classified and named by their shape. Here the surface depends on the signs of and , and the considered coordinates system. Similar analysis for bulk media has been performed to define 2D conical waves. [START_REF] Malaguti | Twodimensional envelope localized waves in the anomalous dispersion regime[END_REF][START_REF] Malaguti | Envelope localized waves of the conical type in linear normally dispersive media[END_REF] Table 1 summarizes the distinct regimes and associated quadrics of rectilinear motion-invariant wavepackets. It reveals that the nature of quadrics is not fully driven by the considered dispersion regime (normal or anomalous).

Quadric surface

Coordinate system

Hyperboloid one sheet Cartesian

Hyperboloid two sheets Cartesian

Elliptic cone Cartesian

Ellipsoid Cartesian

Elliptic paraboloid Cartesian

Paraboloid cylinder Cylindrical

Hyperbolic cylinder Cylindrical

Elliptic cylinder Cylindrical

Table 1: Rectilinear motion-invariant wavepackets in bulk media.

-Case 1: normal dispersion ( ). In the normal dispersion regime, the common situation corresponds to ( close to ). Figure 1(a-b) shows the corresponding surfaces for negative and positive values of in Cartesian coordinates (i.e., in the space), namely hyperboloids of one or two sheets. When this term is null, we find an elliptic cone as shown in Fig. 1(c). As shown in subplots (a-c), the term mainly governs the existence of involved frequencies around the central wavelength of the invariant ST wavepacket. This effect is also clearly visible on the limiting case of ST wave sheets [see conic sections defined as the intersection of the quadric with the plane in Fig. 1(f-h)]. Such programmable correlations between the spatiotemporal degrees of freedom have recently been studied in detail. [START_REF] Kondakci | Diffractionfree space-time light sheets[END_REF][START_REF] Kondakci | Optical spacetime wave packets having arbitrary group velocities in free space[END_REF][START_REF] Bhaduri | Anomalous refraction of optical spacetime wave packets[END_REF][START_REF] Yessenov | Classification of propagation-invariant space-time wave packets in free space: Theory and experiments[END_REF] Such waves take the forms of different hyperbola in the and planes with baseband or passband features, while the relation between and is simply linear, whose slope corresponds to the chosen group velocity. In those subplots, in the vicinity of conic sections, we also point out the local curvature of phasematching induced by dispersion, in strong contrast to standard light cones [START_REF] Kondakci | Diffractionfree space-time light sheets[END_REF][START_REF] Kondakci | Optical spacetime wave packets having arbitrary group velocities in free space[END_REF][START_REF] Bhaduri | Anomalous refraction of optical spacetime wave packets[END_REF][START_REF] Yessenov | Classification of propagation-invariant space-time wave packets in free space: Theory and experiments[END_REF] restricted to free-space propagation of narrowband wavepackets. Next, it is worth mentioning that the sign of parameter can become negative when . As a consequence, one can define a critical group velocity below which , and we then retrieve the typical quadric surface usually obtained in the anomalous dispersion regime (i.e., an ellipsoid) investigated below. Finally, in the case of cylindrical coordinates, note that quadric surfaces [in the space] naturally embed the typical patterns of Xwave (or conical wave) solutions previously studied in bulk [START_REF] Faccio | Conical Waves, Filaments and Nonlinear Filamentation Optics[END_REF][START_REF] Hernandez-Figueroa | Localized Waves[END_REF] or guiding [START_REF] Kibler | Discretized conical waves in multimode optical fibers[END_REF][START_REF] Tarnowski | Numerical modelings of ultrashort pulse propagation and conical emission in multimode optical fibers[END_REF] media, obtained as the intersection of the quadric surfaces and the plane defined by (see Table 1).

-Case 2: anomalous dispersion ( ). In that dispersion regime, the common situation corresponds here to . The associated quadric is an ellipsoid in Cartesian coordinates, as shown in Fig. 1(d) or an elliptic cylinder in cylindrical coordinates. In the latter case, one retrieves the typical signature of O-waves (light bullets) previously studied in bulk media [START_REF] Faccio | Conical Waves, Filaments and Nonlinear Filamentation Optics[END_REF][START_REF] Hernandez-Figueroa | Localized Waves[END_REF] . Again, the limiting case of ST wave sheet can be obtained for in Fig. 1(i). This wave takes the form of distinct ellipses in the and planes, while the linear relation between and is again driven by the chosen group velocity. As previously, the sign of parameter can nevertheless change by strongly varying . Here, one can define a critical group velocity beyond which , and the corresponding quadric becomes a hyperboloid (typical feature of normal dispersion regime above-studied). Such cross-overs from one to another type of surfaces shows that, by carefully choosing , one can easily overcome limitations in the design of ST wavepackets related to the dispersion of the medium under study.

Finally, a specific configuration of light structuring for invariant ST wavepackets is obtained when , again whatever the dispersion regime. Equation (4) becomes linear in terms of frequencies . The associated quadric of rank 2 is an elliptic paraboloid in Cartesian coordinates, as shown in Fig. 1(e) or a paraboloid cylinder in cylindrical coordinates. In the limiting case of ST wave sheet in Fig. 1(j), the conic section degenerates into a simple parabola.

Galilean-transformed (uniform motion and rotation) invariant wavepackets.

One can push further the classification of invariant ST wavepackets by noticing that rectilinear motion-invariant fields form a subgroup of the family which is not only invariant under uniform motion but also invariant at a given rotation around the axis. The latter, called spatiotemporal helicon wavepackets, [START_REF] Faccio | Conical Waves, Filaments and Nonlinear Filamentation Optics[END_REF] belong to the kernel of the space-time screw axis symmetry differential operator , where is an arbitrarily chosen constant and is the z-component of the angular momentum operator. In the following, the study of ST helicon wavepackets is treated in the more general framework of cylindrically symmetric waveguides, i.e., in media whose refractive index may exhibit a radial dependence. The bulk case can then be seen as the limiting case of waveguide of infinite dimension.

For simplicity, we will restrict our study to a scalar approach by considering the weak guidance approximation. The partial differential equation driving the linear propagation of an electric field is then given in cylindrical coordinates by: , (5) where is the radial-dependent refractive index at . By using the separation of transverse and longitudinal variables, the linear evolution of the field along the propagation axis is given in the modal basis [START_REF] Béjot | Spatiotemporal Helicon Wavepackets[END_REF] by: ,

where is the propagation constant of the OAM mode at the frequency . Recall that and refer here to azimuthal and radial indices, where is the topological charge, related to the phase front of OAM modes. The modes form an orthonormal basis so that any electric field can be expressed as , (7) where is the transverse envelope of the mode and are the electric field coordinates in the modal basis. Looking for fields that belong to the kernel of immediately implies that the decomposition of such wavepackets only embeds a family of eigenvectors whose eigenvalues respect:

.

Accordingly, any electric field built from a given family is a diffraction-and dispersion-free ST wavepacket propagating at the group velocity whose intensity continuously rotates around the propagation axis with a spatial period . Before studying how the moderesolved spectrum of such invariant fields is structured in the modal space , one has to emphasize that, generally speaking, the propagation constant in a waveguide (nonlinearly) depends on the topological charge . Hence, still considering only second-order dispersion around a given arbitrarily chosen frequency , the propagation constant of a waveguide can be well approximated as

, ( 9
)
with  is a three-dimensional second-order polynomial:

where the coefficients depend on the geometry of the considered waveguide. For the case of homogeneous bulk media, all coefficients vanish except with being replaced by .

Note that the propagation constant usually does not depend on the sign of . This will have some consequences in the geometric nature of the surfaces describing spatiotemporal helicon wavepackets. Inserting Eq. 9 in Eq. 8, it follows that spatiotemporal helicon wavepackets are contained in surfaces whose algebraic expression are , (10) First, since depends on , it is straightforward to notice that the surface embedding spatiotemporal wavepackets will be composed of two different quadrics (one for and another for ). However, it appears that the rank and nature of the two involved quadrics is the same. The only difference between the two will be the position of their respective principal axis. Then, it is sufficient to study the nature of the quadric for positive only. To this purpose, one has to study the associated matrix to relation ( 10): (11) The relative signs of the eigenvalues of completely classify the nature of the quadric surface in the frequency-resolved modal space.

As an example, we investigate below ST wave structures in a simple and analytical case at the frontier between bulk media and fibers, namely a dispersive medium of finite transversal dimension (radius ). Thus, the propagation constants can be simply described through , where is the root of Bessel function of first kind . Using the formalism described above, with a second-order polynomial fit of (independent of ),

one can find the following matrix :

It appears that analytical formula for eigenvalues of such a matrix of the silica rod remains too complex to be presented here. However, in order to get a simple overview of the nature of quadrics, we considered small values of . This allows to approximatively find the following eigenvalues: .

Accordingly, the quadric surface associated to ST helicon wavepackets in finite bulk media has a rank of 3 (unless ). Moreover, since , the quadric surfaces for and are necessarily either one-or twosheet hyperboloids or elliptic cones in this approximated framework. In the limiting case of (infinite) bulk media, one can easily find tractable analytical formula of eigenvalues from matrix associated to relation (10) with .

In the general framework of matrix (12), typical examples of quadrics are shown in Fig. 2(a-b) in both dispersion regimes. Note that for a clearer illustration, surfaces are plotted even if only discrete values can be found for parameters . We confirm that for each configuration the global surface is composed of two different quadrics (one for and another for ), namely hyperboloids of one or two sheets. Discrete values are then shown for limited cases corresponding to conic sections when or . In the limiting case of infinite bulk media, the quadrics involved for structuring ST helicon wavepackets are theoretically hyperboloids, however, in practice one eigenvalue tends to zero so that surfaces approaches an elliptic/hyperbolic paraboloid, depending on the respective value of , see Fig. 2(c-d). Note that for , we retrieve an elliptic/hyperbolic cylinder as described in Table 1. Numerical construction of wavepackets. Next, we present some examples illustrating our general approach of synthesizing rotating ST wavepackets through quadric phasematching of their spatiotemporal optical components (Eq. 10). In particular, we construct such ST wave structures in the simple case of a fused silica rod with radius µm, both in normal and anomalous dispersion regimes. As a result, for each couple , one must find the roots of the quadric, thus corresponding to frequencies that depict the family of phase-matched modes . Figure 3(a) shows the overall threedimensional pattern of phasematching obtained in the normal dispersion regime, when considering a central wavelength 800 nm, and the arbitrary choice of the constants and (i.e., propagation constants of the wavepacket are those of ). The rotation period is here fixed equal to cm. The linear construction of a helicon wavepacket then results from the superposition of the calculated frequencies . Figure 3(b) presents an example corresponding to a particular case of selected modes with various radial and angular indices , and ), and with equal spectral amplitudes and phases (see circles in Fig. 3a). The successive subplots display the isosurface of the spatiotemporal intensity at half-maximum at distinct propagation distances. The wavepacket looks like a corkscrew in space-time coordinates with extreme localization since all frequencies are in phase. It clearly reveals the spiraling trajectory of the intensity pattern while maintaining its initial distribution over propagation. This helicon wavepacket does not disperse and simply rotates during its propagation around the z-and t-axes because of its inherent invariant nature. The direction of rotation is driven by the sign of . Note that complex engineered 3D-patterns can be obtained by a careful selection of multiple higher-order modes at the phase-matched discrete frequencies.

Another example of quadric phasematching obtained in the anomalous dispersion regime is shown in Fig. 3(c), when considering a central wavelength 1800 nm, and the arbitrary choice of the constants and . The rotation period is now fixed equal to cm. Figure 3(d) presents an example of constructed helicon wavepacket obtained from the superposition of modes with the following radial and angular indices:

, and , and with equal spectral amplitudes and phases (see circles in Fig. 3c). Again we retrieve the typical invariant spiraling trajectory of the intensity pattern. However, we here notice the more complex spatiotemporal pattern made of three consecutive high-intensity sub-pulses with distinct spatial arrangements.

Discussion and conclusion.

The concept of invariant ST wavepackets has emerged as a very promising and powerful tool for controlling the propagation dynamics of pulsed electromagnetic fields. By finely tailoring the modal distribution of the spectral components embedded, one can control in large extent both trajectory and velocity of ST wavepackets, but also the way they spread both in time and space, and this, whatever the dispersion regime. This degree of control let envision extremely exciting applications and new all-optical functionalities. Here, by invoking symmetry arguments, we have shown that invariant ST wavepackets refer to three-dimensional objects whose modal spectral distributions lie on quadrics surfaces.

Figure 3.

Examples of helicon wavepackets constructed from quadric phasematching of spatiotemporal components (from Eq. 10) in a fused silica rod with radius µm. (a,c) Quadric surfaces obtained both in normal and anomalous dispersion regimes. Green dots indicate discretized conic sections for each coordinate. Black circles correspond to the selected modes used for the linear construction of helicon wavepackets depicted in (b,d) through iso-surfaces of their spatiotemporal intensity pattern at half-maximum. The successive subplots obtained every confirm the invariant nature of helicon wavepackets over propagation. The dashed blue line indicates the origin . Projections on planes (shadow plots) are also provided for a clear observation of rotation.

With the extremely rapid progresses made for manipulating ST wavepackets, [START_REF] Cao | Non-spreading Bessel spatiotemporal optical vortices[END_REF][START_REF] Yessenov | Spacetime wave packets localized in all dimensions[END_REF] we hope this work will inspire a wealth of applications, such as recent 3D wavepackets with intriguing vortices and singularities, 17,43- 44 and for which ST invariance feature would require further investigations. Beyond this simple but far-reaching concept of structured light for controlling the linear pulse propagation dynamics, note that the study of quadric waves, the general concept embedding conical waves, will also help in understanding the physics underlying the mechanism of supercontinua generation. For instance, a recent numerical work [START_REF] Béjot | Spatiotemporal Helicon Wavepackets[END_REF] has already reported that the supercontinuum generated by a helicon pulse is naturally structured in a quadric wave that respects the underlying symmetry imposed by the nonlinear propagation of the pump. Finally, the formalism of quadric phasematching is not limited to electromagnetic fields and could find applications in various branches of wave physics, [START_REF] Brunet | Experimental observation of azimuthal shock waves on nonlinear acoustical vortices[END_REF][START_REF] Dominici | Full-Bloch beams and ultrafast Rabi-rotating vortices[END_REF][START_REF] Vieira | Optical control of the topology of laser-plasma accelerators[END_REF][START_REF] Pierce | Coiling free electron matter waves[END_REF][START_REF] Asenjo | Nondiffracting gravitational waves[END_REF] such as acoustics, polaritonics, plasma, gravity, and Bose condensate physics, in order to explore the dynamics of three-dimensional spiraling wavepackets with topological properties.

Conflict of interest.

The authors declare no competing financial interest.

Figure 1 .

 1 Figure 1. Quadrics governing (uniform rectilinear motion) invariant ST wavepacket in bulk fused silica. (a-c) Center wavelength at 800 nm in normal dispersion. (d-e) Center wavelength at 1600 nm in anomalous dispersion. (a) Hyperboloid of one sheet for cm -1 and ps/m. (b) Hyperboloid of two sheets for and ps/m. (c) Elliptic cone for and . (d) Ellipsoid for cm -1 and . (e) Elliptic paraboloid for and ps/m. In

Figure 2 .

 2 Figure 2. Examples of quadric surfaces (from Eq. 10) governing phasematching of space-time helicon wavepackets in a fused silica rod with radius µm (a-b) and bulk fused silica (cd). Subplots (a: Two hyperboloids of one sheet, c: hyperbolic paraboloid) correspond to the normal dispersion regime with central wavelength fixed at 800 nm, and , ps/m, cm. Subplots (b: Two hyperboloids of two sheets, d: elliptic paraboloid) correspond to the anomalous dispersion regime with central wavelength fixed at 1600 nm. , ps/m, and cm (6 cm in subplot d). Blue lines/circles indicate conic sections regarding the limiting cases of non-rotating space-time wavepackets (without superposed OAM) when . Red lines/crosses indicate the simple cases when or .

Funding Sources. The authors acknowledge support from French programs "Investments for the Future" operated by National Research Agency (ISITE-BFC, contract ANR-15-IDEX-03; EIPHI Graduate School, contract ANR-17-EURE-0002; EQUIPEX+ Smartlight, contract ANR-21-ESRE-0040), and from Région Bourgogne Franche-Comté.


C Lacour
email: lacour@onera.fr

Y Maday
email: maday@ann.jussieu.fr

TWO DIFFERENT APPROACHES FOR MATCHING NONCONFORMING GRIDS: THE MORTAR ELEMENT METHOD AND THE FETI METHOD *

Keywords: AMS subject classification: 65F30, 65M60, 65Y05 Domain decomposition, mortar finite element method, nonmatching grids, saddle-point problem, hybrid methods

When using domain decomposition in a finite element framework for the approximation of second order elliptic or parabolic type problems, it has become appealing to tune the mesh of each subdomain to the local behaviour of the solution. The resulting discretization being then nonconforming, different approaches have been advocated to match the admissible discrete functions. We recall here the basics of two of them, the Mortar Element method and the Finite Element Tearing and Interconnecting (FETI) method, and aim at comparing them. The conclusion, both from the theoretical and numerical point of view, is in favor of the mortar element method.

Introduction.

Domain decomposition method has become an important tool for solving large scale problems. One main reason for this is the ease of conception of the algorithms in a framework of parallel implementations. The literature in this field is so large that it is quite impossible to present an overview. The proceedings of the conferences on domain decomposition methods [11, 151 are, maybe, best suited for this.

In addition to the previous reasons it seems interesting to take benefit of the potential of the domain decomposition algorithms to gain even more in flexibility of the discretization. In this direction the ability of using local discretization parameters best fitted to the local characteristic of the solution to be approximated seems natural. The idea of building independently the meshes on each subdomain with no constraint on the interface between the subdomains comes also rapidly into mind. This is particularly appealing when finite element methods are used since it allows then to adapt the mesh locally with no propagation of the refinement from one subdomain to an adjacent one. When the problem to *Received October 1996. Revised February 1997. be solved is of elliptic or parabolic type, in particular, this leads to nonconforming approximations since the meshes are generally nonmatching implying that the discrete spaces cannot be included in the corresponding continuous spaces. The question that has to be addressed is then the gain in CPU (thanks to the non-propagation of refinement in region where it is not necessary) with respect to the actual accuracy (that could be polluted from the nonconformity of the approximation). Another frame of application is related to the treatment of non-constant coefficient since depending on the values of these coefficients, some problems can become stiffer requiring a finer mesh. The nonconforming methods allow again to use the proper discretization parameters best suited to the local characteristic of the problem.

In this paper we recall two basic approaches to deal with nonconformity. These are the Mortar Element method [START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method[END_REF][START_REF] Belgacem | Coupling spectral and finite element discretizations for second order elliptic three dimensional equations[END_REF] and the Finite Element Tearing and Interconnecting (FETI) method [START_REF] Farhat | Using a reduced number of Lagrange multipliers for assembling parallel incomplete field finite element approximations[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. The framework we use for the presentation of the methods is the primal hybrid formulation see [START_REF] Raviart | Primal hybrid finite element methods for second order elliptic equations[END_REF] and [START_REF] Thomas | Formulation mixte des dquations aux ddrivdes partielles du second ordre elliptiques[END_REF]. The parallel between the two methods is given and for both methods, results of numerical analysis and of numerical simulations are given that prove the optimality and the accuracy of the mortar element method.

As is often the case, we shall present the two methods on a simple but generic problem. We shall assume that the domain of computation is a two dimensional polygon, 12 C ~2, and that we want to solve the Poisson problem:

Find u such that

-Au = f, in (1.1) u = 0, over0it
where f is given in L2(~). For both methods many extensions such as to Navier-Stokes equations, e.g., [START_REF] Abdoulaev | The numerical implementation of the domain decomposition method with mortar finite elements for a 3D problem[END_REF][START_REF] Cazabeau | Mdthodes multi-domaines pour la rdsolution des dquations de Navier-Stokes en simulation numdrique directe[END_REF], elasticity problems, e.g., [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Tallec | Three-dimensional domain decomposition methods with nonmatching grids and unstructured coarse solvers[END_REF], have been made that attest that they are not restricted to such a simple problem. The variational statement of problem (1.1) is the starting point of the discretization that consists in rewriting the problem as follows:

Find u E Hi(it ) such that

(1.2)

Vv e Hl(gt), s163

This formulation is very handy for introducing nonoverlapping domain decompositions. Indeed, assume that 12 is partitioned into nonoverlapping (Lipshitz) subdomains This leads naturally to introduce the notation Fk,e = 012 k fq Oft e.

The constraint across the interface Fk,e can be relaxed by inducing the definition of a Lagrange multiplier in the Euler equation. The Lagrange multiplier belongs to a closed subspace M of HI<k<e<KH-1/2(Fk,e). We refer to [START_REF] Thomas | Formulation mixte des dquations aux ddrivdes partielles du second ordre elliptiques[END_REF] and [START_REF] Raviart | Primal hybrid finite element methods for second order elliptic equations[END_REF] for the exact definition of M. Let us simply denote by X* the product of spaces Xk, the problem (1.2) is then equivalent to the following one:

Find u* 9 V* such that K K (1.4)
Vv* E V*, k=lEfd~k VukVvkdX : k~.. :l/~ k and also to the one in primal hybrid formulation over the space X* = IIa= 1K Xk

Find (u*,)~) 9 X* x M such that (1.5) Vv* E X*, V# E M, K = k l<_k<g<K Fk,~ k=l k l(kK~.~_K JFk,l
The Mortar Element method and the FETI method are two methods that allow for designing a discrete version of this problem when different variational approximations are used on each domain f~k. As a result of possibly nonmatching grids, the equality of the traces of the discrete functions across the interfaces has to be relaxed. Starting from the hybrid formulation, it is natural to discretize not only the space X* but also the space of Lagrange multipliers M. Since we face a mixed formulation, a standard requirement is that the space of discrete Lagrange multipliers is not too rich with respect to the space of approximate functions of X*. The two methods that we shall recall differ from the choice of the discrete space for M. REMARK 1.1. There are two other recent methods, at least, that generalize the mortar element method: the three fields formulation introduced in [START_REF] Brezzi | A three fields domain decomposition method[END_REF] and the hybrid approach of [START_REF] Thomas | Une mdthode d'dldments finis hybrides en ddcomposition de domaines[END_REF]. We refer to [3] for a presentation of these methods together with a discussion of these different approaches. REMARK 1.2. In order to have a unified framework, we have chosen to present both methods in a hybrid formulation even though it appears that it is not the best suited for the more general mortar method (e.g., coupling with spectral element method [3]).

Discretization.

We discretize the problem (1.5) by the Galerkin method. Let us consider a parameter h standing for a discretization parameter. For any value of h, for any k, 1 < k < K, we introduce a finite dimensional subspace Yh k of H 1 (~k) M C O (~k) and denote by Xh k the space Yh k M Zk. For any k, 1 < k < K, F k'j, 1 <_ j < j(k) stand for the (eventually curved) segments which coincide with the edges of i2 k, (j(k) denote the number of edges of f~k). We then define the skeleton S as the g ] ]j(k) Fk,J Finally, we choose a union of all edges of all subdomains: S = Uk=l ~j=l finite set A/[ of pairs m = (k, j) such that the F k,j are disjoint from each other. We denote by .ym, and we call mortars, these F k'j.

To describe the discrete space, we begin by defining trace spaces. 9 First, for any k, 1 < k < K and for any j, 1 < j < j(k), we set W~ 'j 'j = y:} 9 Next, for any m* --(k,j) not in A4, we choose a space VVj m of discrete functions on the non-mortar sides, the different choices for the definition of this space lead to different methods in particular the two methods that are presented here: the Mortar Element method and the FETI method.

The product of all these spaces provides a global discretization IIdh of the functions on the skeleton S by

For any m e A/t, we denote by W~ the space W: (m)'j(m). The mortar space is defined by Wh = {~; ~l-~-~ 6 W~,m 6 f14}. The discrete space Vh is the space of functions Vh on f~ such that:

9 For any k, 1 < k < K, Vh, k ----Vhi~k E X k.
9 There exists a function ~ 6 Wh such that: The discretized variational formulation of (1.5) is:

If F k'j is a mortar, Vh,klrk.j ----~. If F k'j is not a mortar
Find u h E Yh such that K (2.1) ~/vh E Yh' k~l k ~Uh'k'~Vh'k K dx=k~l ~= k fVh,kdx.
The problem (2.1) can be reformulated into a saddle point problem. Let ah be the symmetric bilinear form on Xh • Xh:

K ah(Uh,Vh) = ~ f VUh,k.VVh,k dx, J~ k k~-I
and bh the bilinear form on Xh • 17Vh:

(2.2) bh(Vh,#h) = ~ fr (Vh,k --Vh,~),h. l<k<~<K k,t
We can associate to ah the linear operator Ah and to bh the linear operator Bh such that ah(Uh, Vh) -= (AhUh, Vh) and bh(Vh,Ph) = (BhVh,#h). Therefore, the problem (2.1) admits a following saddle-point formulation:

Find the pair (Uh,)~h) in Xh X 17Vh such that Ahuh + B~Ah = fh

(2.3) BhUh ----O.
3 The mortar method.

The mortar method corresponds to a particular choice for the gluing functions r over the non-mortar sides. We define each I~ as a proper subspace of codimension 2 in W k'j. In case of first order approximation, it is the subspace of all functions with zero slope on the elements including the end points. It is then an easy matter to make the following choice of basis DEFINITION 3.1. In case of first order finite element spaces Xkh the basis functions -h of CV k'j are chosen as piecewise linear functions with zero slopes at the end points of interfaces, that vanish at any inner point of the discretization ofrk,J except one (see Fig.

The interface matrix.

In this section, we will explain how to build the interface matrix Bh which makes the correspondence between the degrees of freedom on interface and the degrees of freedom of the Lagrange multipliers. In all what follows, we consider the case of piecewise linear finite elements. We focus our attention on an interface ~/= Xox]y = ~,y --/3[ between two subdomains ~t 1 and ~2 and we assume that the side of ~2 is the mortar one. We have to glue two subdomains with nonmatching grids on their common interface. Therefore, we have to solve:

ff (ul - = 0 1
which can equivalently be written as

U 1 = O. 1 B1 ~32
We consider the shape functions of ui piecewise linear and the shape functions of r as already introduced in definition 3.1. Since the non-mortar side corresponds to domain 1, the gluing function r is defined on the grid of ~t 1

Computing B1.

We denote by n the number of vertices of the triangulation of fll that belong to the interface. Basis functions corresponding to r are hk with k from 2 to n -1 and basis functions corresponding to ul are hi. 9 Ifk=2=:~h2=hl+h2. 9 Ifk=3,...,n-2~-hk=hk. 9 Ifk=n-l~hk=hn-l+hn. The evaluation of B1 consists in computing:

ia,, ~ ul(ai)hi-hk

2 to n-1 .
,& with k from

1 i=1
We remark that B1 is a rectangular matrix with n -2 lines and n columns, so that the inner values ul(ai)i=2 ....... 2 are entirely determined from the values of uz(al), ul(an) and the complete set u2(y), y E (A1, ..., Am).

How to implement BI?

9 End points.

1. For i---1 0 ifj > 3.

= 0 ifj>2, = hi = __a2-al = __a2-al 1
2. Similarly for i = n -1.

9 Inner points.

~ h&

fo. ai + l hi-hi Basis functions corresponding to r are the same as before hk for k = 2,..., n-1 and prolongated by 0 over F 2,m for some m in case F 1,n does not coincide with F 2,m. We set F 2,m = ]a, ~I-Basis functions corresponding to u2 are Hi as is illustrated in Fig. 3.2.

a 1 i- = 0 ifjr i, i+1 3
The technique now consists in introducing a finer grid of points (~/)t=l ..... F on the interface as shown in Fig. 3 As for the mortar approach for the Lagrange multipliers, we have to glue two subdomains with non-matching grids on their interface. So, we have to satisfy the condition ~a an Computation: The integral can equivalently be written as:

r - u2r = 0. 1 y 9 Y

B1 B2

The evaluation of each term is decomposed over each element ]ai, ai+l[ of the interface l(zo_ r

gl = ~ + i=2 ~-1 ,i a i (4.3) = ~(/_111-Fx Ix+l, ] ai-ai_ldx i=2 ~ r L---~--(ai -ai-1) + ai-1 2 + f_: 1- x_V r [--V -ai) +ailai+12-aidx )--
and the integrals over each element are computed exactly by using a Gauss-type quadrature formula. A similar treatment is done for •2 based on the points Ak.

The basis functions that are used to span the space of polynomials are the set of Legendre polynomials Ln.

5 Numerical experiments.

Generalities.

For mixed formulations, as the ones that are encountered in this type of method, it is well known that a crucial argument is related to the compatibility between the discrete space of Lagrange multipliers and the space of discrete subdomain solutions. In order to have a unique solution, it is important that the space of Lagrange multipliers is not too rich. This is expressed by the standard inf-sup condition recalled in the following hypothesis. where the constant/3~ is > 0 but may depend on the discretization parameters.

Of course, in order to have an optimal approximation, this dependency of the inf-sup constant has to be as low as possible, the best being that /3~ is independent of the discretization parameter. For the Mortar method, it is proven that it is the case: the inf-sup condition is independent of the discretization parameter [START_REF] Belgacem | The mortar finite element method with Lagrange multipliers[END_REF] and the optimality of the discrete method results. Actually, it is proven in [START_REF] Belgacem | The mortar finite element method with Lagrange multipliers[END_REF] that the inf-sup condition holds on a smaller space for the v. Let us denote by Zh the subspace of Xh of elements that vanish on the mortars (and not on the nonmortars), then the following inf-sup condition is satisfied

(5.1) infsup ~-~(k'~)fr~'~vkh#hds _>/3~ ((k~,~) :_~_ )1/2
REMARK 5.1. Any choice of space for the Lagrange multiplier as being a subspace of the original proposed space results also in a constant inf-sup condition. This is nevertheless at the price of a degradation of the approximation of the Lagrange multiplier s that actually pollutes the global solution.

For the FETI method, no theory is yet available but it is noticeable from the numerical point of view that the inf-sup constant is not independant of the discretization parameter, i.e the degree of the local Lagrange multipliers as is presented in the next section.

Numerical evaluation of the inf-sup condition.

This part deals with the relation between the eigenvalues of the algebraic system resulting from the implementation of the Shur dual method in the case of nonconforming grids and the inf-sup constant arising in the verification of the inf-sup condition. As already noticed in [START_REF] Maday | Analysis of iterative methods for the steady and unsteady Stokes problem: Application to spectral element discretizations[END_REF], the elements v~ E Xh that realize the supremum in this inf-sup condition, i.

Let us first consider the mortar element method. Let ]/h be given in

VPe~ -II,uh(p)ll_}

where #h (P) denotes the element of Wh the nodal values of which, by definition, are equal to P. Recalling that the matrix BtkA-~lBk is symmetric, it follows that inf

E(BtkA;1BkQ, Q) : O/mi n QC~tM'IIQIIe~(~M)=I k
where amin is the smallest eigenvalue of t -1

~'~k BkAk Bk. It is well known also from the equivalence of the norms over Wh given by the condensation of the mass matrix that, for any vector Q in ~M,

"Q"e2(~M) ~--i~k, hk"Ph(O)l'2L2(rk,~ )
moreover, the standard imbeddings and the inverse inequalities over the finite element functions give REMARK 5.2. Actually it is reasonable to think that for elements that realize the infimum in the inf-sup constant, the inverse inequality (5.2) scales like h -1/2, meaning that (in case, for simplification, where hrnin and hm~x are comparable) fq m__in (5.3) t3~ ~_,JV h-2 "

(5.2) [I,h(Q)ll(Ho 89 < II,h(Q)llL~(r~.~) _< h;U211,h(Q)ll(Hoto(r~.m .
It is also reasonable to think that C~ma x scales exactly as a constant times h (indeed it is lower bounded by the value of frk.~ Vh#h ds/(]lVh [[.[[phl[) obtained for ~h being an approximation of a sin wave and Vh being an approximation of an harmonic extension of the sin wave).

These two heuristic statements come in accordance with the numerical experiments reported in Figs. 5.3, 5.4, where amin is proven to scale like h 2 and ~ as c/h.

Let us now consider the polynomial FETI approach. If the Legendre basis is chosen to represent the Lagrange multiplier, the only difference with respect to the previous analysis is in the relation between the H-1/2-norm of #N and the ~2-norm of its coefficients. Since this g2-norm is equivalent to the L2-norm of #g, the only element comes from the inverse inequality (5.2) yielding g 2

_<c--

The numerical experiment reported in Fig. 5.2 shows that the condition number a grows exponentially fast with the degree N of the polynomial Lagrange multiplier so that the inf-sup constant has a very bad behaviour. 3. The influence of this bad inf-sup condition is not only noticeable on the approximation result but yields also a degradation of the performance of the algorithms for solving the problem. In order to circumvent this drawback it has been proposed in [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] to use "piecewise (relatively) low order polynomial" instead of full polynomials. The FETI method then gets closer to the mortar method.

In addition, the verification of the inf-sup condition is not certified and is a prerequisite of the method. To verify this condition, the Lagrange multiplier space should not be too rich. But, the less rich the space is, the less accurately the relation Biui -B2u2 = 0 is verified. In practice, it is difficult to check whether this condition is satisfied or not. Therefore, it is recommended to start with low order approximations of ~ and refine them. The rule of thumb consists in increasing the degree of polynomial N until N -1 is less than the number of inner degrees of freedom.

We compare now the Mortar Element method with the FETI method with polynomial Lagrange multipliers. We study the problem given in (1.1) on a square with f = 2 sin(x) sin(y). The analytical solution is then u = sin(x) sin(y). In Fig. 5.1 we have plotted a section of the numerical solution obtained with the same finite element methods over each subdomain and with different Lagrange functions. The domain is a square, decomposed into 4 (2 x 2) squares. Each square is provided with a quadrangulation: the two upper with 24 x 24 elements, the two lower with 32 x 32 elements. The mortar results are obtained by using the choice advocated in section 3 (the gluing functions are on the coarser mesh, i.e., with 24 elements) and the polynomial results are obtained by using local polynomials with degree 23 so that the dimension of the Lagrange space is the same. (As an illustration of the rule of thumb, we increase the degree of polynomial until N = 23 and for N = 24 we remark that the system is not invertible.)

We do not pretend that this simple numerical example is a definitive statement on the compared qualities of the FETI and the Mortar approaches. Nevertheless, it is in concordance with the evaluation of the inf-sup condition. What we claim overall is that the numerical results presented in this section constitute a definitive answer as regards the evaluation of the inf-sup condition and the bad condition number of the resulting algebraic system.

6 Conclusions.

In this paper, we have compared two different approaches for the matching of nonconforming finite element methods: the Mortar method and the polynomial FETI method. We have presented them in a unified version and focussed on the differences that reduce to be only on the choice of the Lagrange multiplier space that allow to glue together the different meshes. We have pointed out that the choice of mortar method is more appropriate as it allows to 9 certify that the compatibility condition between the discrete spaces is satisfied 9 provide an inf-sup condition that is independent of the discretization parameter 9 lead to an algebraic system with well-conditioned matrices. This is the reason why we have preferred this choice for the generalization to more complex situations. The already good condition number of the matrix is improved by developing preconditioners in the same spirit as the one proposed in [START_REF] Roux | M~thode de dgcomposition de domaine it l'aide de multiplicateurs de Lagrange et application it la rdsolution en parall~le des dquations de l'dlasticitd lindaire[END_REF]. Note that a hierarchical basis for the mortar space has also been introduced that improved even more the convergence of the iterative resolution of the problem. We refer to [START_REF] Yserentant | On the multi-level splitting of finite element spaces[END_REF] and [START_REF] Lacour | Iterative substructuring preconditioners for the mortar finite element method[END_REF] for the presentation of these techniques.

  .2) can be rewritten as follows:Find u E H] (it) such thatK K fa = ~ flakvlak dx Vv e Hi(it), E V(Ul~k)V(Vl ~k)dx Z k k=l k k=l Instead of searching an element u defined globally over ft, it is more convenient, especially when local discretizations are to be used, to search for a K-uple u* = (Ul,..., ug). The space V* spanned by these restrictions V* = {v* = (vl,...,vg), 3v 9 H~(ft),Yk, 1 < k < K, Vk = Vlnk} can be conveniently rewritten as an aggregate of the local spaces Xk = {vk E Hl(ftk), vk = O over Of~ k fq cof~} as follows V*={v*=(vl,...,vg) 9 Vk, e, l < k,e <_ K, vk = v~ over Of~ k N Ofte}.

2 Figure 2 . 1 :

 221 Figure 2.1: Nodal basis of Lagrange multipliers space; the points ai are the vertices of the triangulation of ~k on F k'j

  (2.1)).

Figure 3 . 1 :

 31 Figure 3.1: Lagrange multipliers space.

1 2 2 2. 1 n

 21 Similarly, for i = n il ~ h~hj : 0 ifj<n-1, = hn --an -an-

3 3 ~al

 3 ~-: hi-hi+l-~ /_11 (~--~)(~"~-)ai+12--ai

  .

4 .~Figure 3 . 2 :Figure 3 . 3 :Figure 3 . 4 :

 4323334 Figure 3.2: Interface matrix calculus.

Figure 4 . 1 :

 41 Figure 4.1: Polynomial Lagrange multipliers.

HYPOTHESIS 5 . 1 (

 51 LADYZHENSKAYA, BABUSKA AND BREZZI). Let us denote by Fk,~ any non-mortar edge of the subdomain ~k, we have:

  to deduce from (5.1) that another inf-sup condition then holds inf sup }-~(k,e) Jr~,e h --.hew,, .,,,~x,, IIv,,ll.->/~ I1~"11~ 89 '(Hoo(Pk,e))

V

  the element ~ of Xa solution of the problem s f k k (k,~) JFk,e where the solutions are locally chosen to be zero average over the flotting subdomains in order to ensure the uniqueness of the solution ~. It follows that sup ' = (~k _ h Jt~h ds (k,t) Fk,~ If P denotes the vector in ]~M of nodal values of #h in the basis h, then it is immediate, from the previous definitions to note that E (BtkA-~' BkP' P) E f ((:;k -(k,e) (Hoo (Fk'e)) This provides the relation r (v k -V~h)#h ds E(k,f) JFk,~ h inf sup .h~Wh v,,~x,, Ilvhll,ll,~hll-89 = ,/ inf Y~'k(BtkAklBkP' P)

Finally, summing upFigure 5 . 1 :

 51 Figure 5.1: Solution with the Mortar method and the FETI method.

Figure 5 . 2 :Figure 5 . 3 :

 5253 Figure 5.2: Behaviour of the Dual matrix' condition number versus the degree of polynomial, same mesh as in Fig. 5.1, the horizontal line corresponds to the mortar Lagrange multipliers.

Figure 5 . 4 :

 54 Figure 5.4: Mortar Lagrange multipliers: behaviour of f~ versus the mesh size, same mesh as in figure 5.1

Behaviour of beta versus the mesh size

		0.1
		0.09
		0.08
		0.07
		0.06
	i	0.05
		0.04
		0.03
		0.02
		0.01
		0 1 H	0,8	0.6	H/2	0.4	H/4	0.2	H/8
					H

Acknowledgement.

The authors want to express their thanks to C. Farhat and F. X. Roux for advice and helpful discussions.

where k goes from 2 to n -1 and where 9 n is the number of points of the triangulation of [21 and the interface. 9 m is the number of points of the triangulation of ~2 and the interface. 9 :f is the number of points on the fine grid.

The integration over ] ~f~f+l[ is approximated by a trapezoidal rule, we obtain:

with k from 2 to n -1. 1

How to implement B2?

and similar relations in other cases.

4 The FETI method.

The FETI method differs from the previous mortar method from the choice of the Lagrange multipliers. Here the skeleton S is further decomposed since we consider the partition S -~-(,.Jl~_k<e<N~k~, where we recall that Fk~ is the segment such that Fke = ~k n-~. Over each entity Fke a space of local polynomials is introduced with degree _< nk,e.


Maxence Grand
email: maxence.grand@univ-grenoble-alpes.fr

Damien Pellier
email: damien.pellier@univ-grenoble-alpes.fr

Humbert Fiorino
email: humbert.fiorino@univ-grenoble-alpes.fr

An Accurate HDDL Domain Learning Algorithm from Partial and Noisy Observations

The Hierarchical Task Network (HTN) formalism is very expressive and used to express a wide variety of planning problems. In contrast to the classical STRIPS formalism in which only the action model needs to be specified, the HTN formalism requires to specify, in addition, the tasks of the problem and their decomposition into subtasks, called HTN methods. For this reason, hand-encoding HTN problems is considered more difficult and more error-prone by experts than classical planning problem. To tackle this problem, we propose a new approach (HierAMLSI) based on grammar induction to acquire HTN planning domain knowledge, by learning action models and HTN methods with their preconditions. Unlike other approaches, HierAMLSI is able to learn both actions and methods with noisy and partial inputs observation with a high level or accuracy.

Introduction

The Hierarchical Task Network (HTN) formalism [START_REF] Erol | HTN Planning: Complexity and Expressivity[END_REF] is very expressive and used to express a wide variety of planning problems. This formalism allows planners to exploit domain knowledge to solve problems more efficiently [START_REF] Nau | Applications of SHOP and SHOP2[END_REF] when planning problems can be naturally decomposed hierarchically in terms of tasks and task decompositions. The standard language used to model HTN problem is HDDL (Hierachical Domain Description Language) [START_REF] Höller | HDDL: An Extension to PDDL for Expressing Hierarchical Planning Problems[END_REF]. In contrast to the classical PDDL language used to model STRIPS problems in which only the action model needs to be specified, HDDL requires to specify the task model of the problem. A task model can be primitive and compound. A primitive task model is described by PDDL operators. A compound tasks model is described using HTN methods. An HTN method describes the set of primitive and/or compound task required to decompose a specific compound task. For this reason, hand-encoding HTN problems is considered more difficult and more error-prone by experts than classical planning problem. This makes it all the more necessary to develop techniques to learn HTN domains.

Many machine learning approaches have been proposed to facilitate the acquisition of PDDL domain acquisition and to learn the underlying action model, e.g, ARMS [START_REF] Yang | Learning action models from plan examples using weighted MAX-SAT[END_REF], FAMA [START_REF] Aineto | Learning action models with minimal observability[END_REF], LOCM [START_REF] Cresswell | Acquiring planning domain models using LOCM[END_REF], LSONIO [START_REF] Mourão | Learning STRIPS Operators from Noisy and Incomplete Observations[END_REF], AMLSI (Grand, Fiorino, and Pellier 2020a,b). In these approaches, training data are either (possibly noisy and partial) intermediate states and plans previously generated by a planner, or randomly generated action sequences (i.e. random walks). On the other hand, few approaches have been proposed to learn HTN domains. However, it is possible to mention CAMEL [START_REF] Ilghami | CaMeL: Learning Method Preconditions for HTN Planning[END_REF], HTN-Maker [START_REF] Hogg | HTN-MAKER: Learning HTNs with Minimal Additional Knowledge Engineering Required[END_REF][START_REF] Hogg | Learning hierarchical task networks for nondeterministic planning domains[END_REF], LHTNDT (Nargesian and Ghassem-Sani 2008) or HTN-Learner [START_REF] Zhuo | Learning HTN Method Preconditions and Action Models from Partial Observations[END_REF]. The major drawbacks of these approaches are: (1) they consider to have complete and noiseless observations as input; (2) they only learn HTN methods except HTN-Learner, i.e., they consider that the action model is known a priori and (3) the learned domains are not accurate enough to be used "as is" in a planner. A step of expert proofreading is still necessary to correct them. Even small syntactical errors can make sometime the learned domains useless for planning

To deal with these drawbacks, we propose in this paper, a new learning algorithm for HDDL domains, called Hi-erAMLSI. HierAMLSI is based on AMLSI (Grand, Fiorino, and Pellier 2020a,b), a PDDL domain learner based on grammar induction. HierAMLSI takes as input a set of partial and noisy observations and learns a full HDDL planning domain with action model and HTN methods. We show experimentally that HierAMLSI is highly accurate even with highly partial and noisy learning datasets minimising HDDL domain proofreading by experts. In many HDDL ICP benchmarks HierAMLSI does not require any correction of the learned domains at all.

The rest of the paper is organized as follows. In section 2 we present the problem statement. In section 3 we give some backgrounds on the AMLSI approach. In section 4, we detail the HierAMLSI steps. Then, section 5 evaluates the performance of HierAMLSI on IPC benchmarks. Finally, Section 6 concludes with the related works.

Formal Framework

Section 2.1 introduces a formalization of STRIPS planning domain learning consisting in learning a transition function of a grounded planning domain and in expressing it as PDDL operators and Section 2.2 extends this formalization to HTN domains.

STRIPS Planning

In this section we use definitions and notations proposed by [START_REF] Höller | Assessing the Expressivity of Planning Formalisms through the Comparison to Formal Languages[END_REF]) and adapt them to the learning problem.

A STRIPS planning problem is a tuple P = (L, A, S, s 0 , s, δ, λ), where L is a set of logical propositions describing the world states, S is a set of state labels, s 0 ∈ S is the label of the initial state, and g ⊆ S is the set of goal labels. λ is an observation function λ : S → 2 L that assigns to each state label the set of logical propositions true in that state. A is a set of action labels.Action preconditions, positive and negative effects are given by the functions prec, add and del that are included in δ = (prec, add, del). prec is defined as prec : A → 2 L . The functions add and del are defined in the same way. Without loss of generality, we chose this unusual formal framework inspired by [START_REF] Höller | Assessing the Expressivity of Planning Formalisms through the Comparison to Formal Languages[END_REF]) in order to define the STRIPS learning problem as the lifting of a state transition system into a propositional language.

The function τ : A × S → {true, f alse} returns whether an action is applicable to a state, i.e. τ (a, s) ⇔ prec(a) ⊆ λ(s). Whenever action a is applicable in state s i , the state transition function γ : A × S → S returns the resulting state

s i+1 = γ(s i , a) such that λ(s i+1 ) = [λ(s i ) \ del(a)] ∪ add(a).
A sequence (a 0 a 1 . . . a n ) of actions is applicable to a state s 0 when each action a i with 0 ≤ i ≤ n is applicable to the state s i . Given an applicable sequence

(a 0 a 1 . . . a n ) in state s 0 , γ(s 0 , (a 0 a 1 . . . a n )) = γ(γ(s 0 , a 0 ), (a 1 . . . a n )) = s n+1 .
It is important to note that this recursive definition of γ entails the generation of a sequence of states (s 0 s 1 . . . s n+1 ). A goal state is a state s such that g ∈ G and λ(g) ⊆ λ(s). s satisfies g, i.e. s |= g if and only if s is a goal state. An action sequence is a solution plan to a planning problem P if and only if it is applicable to s 0 and entails a goal state.

In formal languages, a set of rules is given that describe the structure of valid words and the language is the set of these words. For STRIPS planning problem P = (L, A, S, s 0 , G, δ, λ), this language is defined as (0 ≤ i ≤ n):

L(P ) = {ω = (a 0 a 1 . . . a n )|a i ∈ A, γ(s 0 , ω) |= g}
We know that the set of languages generated by STRIPS planning problems are regular languages [START_REF] Höller | Assessing the Expressivity of Planning Formalisms through the Comparison to Formal Languages[END_REF]. In other words, a STRIPS planning problem P = (L, A, S, s 0 , G, δ, λ) generates a language L(P ) that is equivalent to a Deterministic Finite Automaton (DFA) Σ = (S, A, γ). S and A are respectively the nodes and the arcs of the DFA, and γ is the transition function.

For any arc a ∈ A, we call pre-set of a the set preset(a) = {s ∈ S | γ(s, a) = s ′ } and post-set of a the set postset(a) = {s ′ ∈ S | γ(s, a) = s ′ } (see Figure 1).

A STRIPS learning problem is as follow: given a set of observations Ω ⊆ L(P ), is it possible to learn the DFA Σ, and then infer P ?

For instance, suppose Ω = {a, ab, ba, bab, abb, . . . (Grand, Fiorino, and Pellier 2020b) show that it is possible to learn Σ (see Figure 1)

} such that s 0 a - → s 2 , s 0 a - → s 2 b - → s 2 , s 0 b - → s 1 a - → s 2 , s 0 b - → s 1 a - → s 2 b - → s 2 , s 0 a - → s 2 b - → s 2 b - → s 2 . . .
s 0 start s 2 s 1 b a a b postset(a)
preset(a) Figure 1: An example of DFA with pre-states and post-states and infer P with actions {a, b}, the initial state s 0 and some states marked as goal G = {s 2 }.

HTN Planning

By extension based on the notation of [START_REF] Höller | Translating totally ordered HTN planning problems to classical planning problems using regular approximation of context-free languages[END_REF], an HTN planning problem is a tuple P = (L, C, A, S, M, s 0 , w I , g, δ, τ, λ, σ, ζ). As for STRIPS problems, L is a set of logical propositions describing the world states, S is a set of state labels, s 0 ∈ S is the label of the initial state, g ⊆ S is the set of goal label, λ is the observation function and preconditions, positive and negative effects are given by the functions prec, add and del included in δ.

A is the set action (or primitive task) labels and C is a set of compound (or non primitive) task labels, with C ∩ A = ∅. Tasks are maintained in task networks. A task network is a sequence of tasks (for simplicity, we consider only Totally Ordered domain). Let T = C ∪ A. A task network is an element out of T * ( * is the Kleene operator). Compound tasks are decomposed using methods. The set M contains all method labels. Methods are defined by the function σ : M → C × T * . Then, a coumpound task c is decomposable in a state s if and only if there exists a revelant method m ∈ M such that: σ(m) = (c, ϕ) and prec(m) ∈ s. The function ζ : T * × S → T * gives the decomposition function. For a totally orderer task network ω = ω 1 tω 2 , ζ is defined as follows:

ζ(ω 1 tω 2 , s) =      ω 1 tω 2 if t is a primitive task ω 1 ϕω 2 if t is a coumpound task and t is decomposable in γ(ω 1 , s) ∅ Otherwise
As ω 1 tω 2 is a totally ordered task network, ω 1 contains only primitive tasks. We denote ω → * ω ′ that ω can be decomposed into ω ′ by 0 or more method applications. Finally, ω I is the initial task network.

A solution to an HTN planning problem is a task network ω with:

1. ω I → * ω, i.e. it can be reached by decomposing ω I . 2. ω ∈ A * , i.e. all tasks are primitive. 3. γ(s 0 , ω) |= g, i.e. ω is applicable in s o and results in a goal state. Finally, we can define an HTN planning problem P = (L, C, A, S, M, s 0 , w I , g, δ, τ, λ, σ, ζ) as a formal language:

L(P ) = {ω = (t 0 t 1 . . . t n )|t i ∈ A, γ(s 0 , ω) |= g, ω I → * ω}
An HTN learning problem is as follow: given a set of observations Ω ⊆ L(P ), is it possible to learn P ?

Unlike STRIPS planning problems, the language L(P ) is not necessary regular [START_REF] Höller | Language Classification of Hierarchical Planning Problems[END_REF] and cannot be represented as a DFA. As mentioned by [START_REF] Höller | Language Classification of Hierarchical Planning Problems[END_REF][START_REF] Höller | Translating totally ordered HTN planning problems to classical planning problems using regular approximation of context-free languages[END_REF]), L(P ) is the intersection of two languages:

1. L H (P ) = {ω ∈ A * |w I → * ω} which is defined by the decomposition hierarchy, i.e. by the compound tasks and methods. 2. L C (P ) = {ω ∈ A * |γ(s 0 , ω) ∈ g} which is defined by the state transition system defined by the preconditions and effects of the primitive tasks. This language is regular.

The key idea of our approach is to learn the DFA Σ = (S, A, γ) corresponding to the regular language L C (P ), and modify the DFA to approximate the language L(P ) with the DFA Σ = (S, T, γ) and then infer P .

The AMLSI Approach

In this section we will summarized the AMLSI approach on which HierAMLSI is based. For more detail see (Grand, Fiorino, and Pellier 2020b,a).

AMLSI generates the set of observations Ω by using random walks to learn Σ = (S, A, γ) and deduce P = (L, A, S, s 0 , G, δ, λ). AMLSI assumes L, A, S, s 0 known and the observation function λ possibly partial and noisy (a partial observation is a state where some propositions are missing and a noisy observation is a state where the truth value of a proposition is erroneous). No knowledge of the goal states G is required. Once Σ is learned, AMLSI has to deduce δ from the transition function γ. Concretely, δ can be represented as a STRIPS planning domain containing all the actions of the problem P and by induction the classical PDDL operators.

The AMLSI algorithm consists of 4 steps: (1) generation of the observations, (2) learning the DFA corresponding to the observations, (3) induction of the PDDL operators from the learned DFA; (4) finally, refinement of these operators to deal with noisy and partial state observations:

Step 1: AMLSI generates a random walk by applying an action from the initial state of the problem. If the action is applicable in the current state the sequence of actions from the initial state is valid and is added to I + , the set of positive samples. Otherwise the random walk is stopped and the sequence is added to I -, the set negative samples.

Step 2: To learn the DFA Σ = (S, A, γ) AMLSI uses a variant (Grand, Fiorino, and Pellier 2020b) of a classical regular grammar learning algorithm called RPNI [START_REF] Oncina | Inferring regular languages in polynomial update time[END_REF]. The learning is based on both I + and I -.

Step 3: AMLSI begins by inducing the preconditions and effects of the actions. For the preconditions prec(a) of action a, AMLSI computes the logical propositions that are in all the states preceding a in Σ:

prec(a) = ∩ s∈preset(a) λ(s)
For the positive effects add(a) of action a, AMLSI computes the logical propositions that are never in states before the execution of a, and always present after a execution:

add(a) = ∩ s∈preset(a) λ(s) \ prec(a) Dually, del(a) = prec(a) \ ∩ s∈postset(a) λ(s)
Once preconditions and effects are induced, actions are lifted to PDDL operators based on OI-subsumption (subsumption under Object Identity) [START_REF] Esposito | Multistrategy Theory Revision: Induction and Abduction in INTHELEX[END_REF]: first of all, constant symbols in preconditions and effects are substituted by variable symbols. Then, the less general preconditions and effects, i.e. preconditions and effects encoding as many propositions as possible, are computed as intersection sets. This generalization method allows to ensure that all the necessary preconditions, i.e. the preconditions allowing to differentiate the states where actions are applicable from states where they are not, to be rightfully coded in the corresponding operators.

Step 4: To deal with noisy and partial state observations, AMLSI starts by refining the operator effects to ensure that the generated operators allow to regenerate the induced DFA. To that end, AMLSI adds all the effects ensuring that each transition in the DFA are feasible. Then, AMLSI refines the preconditions of the operators. As in [START_REF] Yang | Learning action models from plan examples using weighted MAX-SAT[END_REF], it makes the following assumptions: the negative effects of an operator must be used in its preconditions. Thus, for each negative effect of an operator, AMLSI adds the corresponding proposition in the preconditions. Since effect refinements depend on preconditions and precondition refinements depend on effects, AMLSI repeats these two refinements steps until convergence, i.e., no more precondition or effect is added. Finally, AMLSI performs a Tabu Search to improve the PDDL operators independently of the induced DFA, on which operator generation is based. Once the Tabu Search reaches a local optimum, AMLSI repeats all the three refinement steps until convergence.

The HierAMLSI approach

HierAMLSI generates the set of observations Ω by using random walks to learn Σ = (S, T, γ) and deduce P = (L, C, A, S, M, s 0 , w I , g, δ, τ, λ, σ, ζ). HierAMLSI makes the same assumptions than AMLSI: L, A, S, s 0 are known but also C, the decomposition function ζ and the observation function λ possibly partial and noisy (a partial observation is a state where some propositions are missing and a noisy observation is a state where the truth value of a proposition is erroneous). No knowledge of the goal states g is required. Once Σ is learned, HierAMLSI has to deduce δ from the transition function γ and σ from the decomposition function ζ.

Like AMLSI, HierAMLSI uses random walks to generate Ω and makes the same assumptions.

HierAMLSI has been devised to solve HTN learning problems. In practice, it computes P as HDDL domain and problem files [START_REF] Höller | HDDL: An Extension to PDDL for Expressing Hierarchical Planning Problems[END_REF][START_REF] Höller | Hierarchical Planning in the IPC[END_REF].

Regarding the training dataset, HierAMLSI uses random walks to generate Ω. HierAMLSI makes the same assumptions than AMLSI. Once the DFA is learned, HierAMLSI generates the set of methods σ and infers the action precondition, positive and negative effect functions in δ from the state transition function γ. Finally, the methods and operators of the HDDL domain file are induced from σ and δ.

The HierAMLSI approach consists of 4 steps: 1. Generation of the observations. HierAMLSI produces a set of observations Ω by using a random walk. In section 4.1, we will present how HierAMLSI is able to efficiently exploit these observations by taking into account not only the fact that some task are decomposable in certain states and their decomposition but also that others are not. 2. DFA Learning. HierAMLSI learns a DFA approximating the language L(P ) (see Section 4.2). 3. HTN Methods generation. HierAMLSI generates from the DFA learned previously a set of HTN Methods allowing to decompose all tasks observed in Ω (see Section 4.3). 4. Action model and HTN Methods precondition learning. Once HTN Methods have been learned, HierAMLSI have to learn the action model, i.e. primitive tasks preconditions and effects and the HTN Methods preconditions. To do this, HierAMLSI treat HTN Methods as primitive tasks and learn an action model containing both methods and primitive tasks using the learning and refinement techniques proposed by the AMLSI approach (see Section 3). The rest of this Section will be illustrated using the IPC 1 Gripper domain. In this domain, a robot moves balls in different rooms using its two grippers r and l. This domain contains three compound tasks: goto(r a ) the robot goes into the room r a , move1ball(b 1 r a ) the robots move the ball b 1 in the room r a and move2balls(b 1 b 2 r a ) the robot moves balls b 1 and b 2 in the room r a .

Observation Generation

The data generation process is similar to the generation method of the AMLSI algorithm (Grand, Fiorino, and Pellier 2020b). To generate the observations in Ω, HierAMLSI uses random walks by querying a blackbox. HierAMLSI chooses randomly a (primitive or compound) task t. If the task t is decomposable, the blackbox returns the final decomposition ϕ containing only primitive task and this decomposition is added to the current primitive task sequence. This procedure is repeated until the selected task is not applicable to the current state. The applicable prefix of the primitive task sequence is then added to I + , the set of positive samples, and the complete sequence, whose last task is not applicable, is added to I -, the set of negative samples. Random walks are repeated until I + and I -achieve an arbitrary size.

DFA Learning

As mentioned in Section 2 the language L(P ) is not necessary regular, then the purpose of this step is to learn a DFA 1 https://www.icaps-conference.org/competitions/ Algorithm 1: Heuristic HTN Methods learning

1 Initialization(M ) 2 for i = 1 : |C| do 3 for c ∈ C do 4 M ′ [c] ← greedy(c, M, i) 5 for c ∈ C do 6 if |M ′ [c]| < |M [c]| then 7 M [c] ← M ′ [c] 8 return M
approximating this language. More precisely, the DFA learning step is divided in 2 steps: (1) HierAMLSI learns a DFA corresponding to the language L C (P ) which is defined by the state transition system defined by the preconditions and effects of the primitive tasks and (2) HierAMLSI adds transitions to represent compound tasks in the DFA to allow to approximate the language L(P ).

Step 1: Primitive task DFA Learning HierAMLSI starts by using the DFA Learning algorithm proposed by (Grand, Fiorino, and Pellier 2020b) to learn the DFA containing only primitive tasks.

Step 2: Task DFA Induction Once the Primitive Task DFA has been learned, HierAMLSI induces the Task DFA by adding compound task transition in the DFA, i.e. by adding transitions whose labels are compound task labels. For instance, suppose we have the compound task

HTN Methods Learning

Once the Task DFA is induced HierAMLSI can directly extract HTN Methods from the Task DFA. However, it is possible that a large number of HTN Methods has been generated. For instance, let's take the Task DFA in Figure 2. For the compound task move1ball(?b ?r) HierAMLSI can generate several methods: ϕ 1 = (move(?r ?r 1 ), pick(?b ?g ?r 1 ), move(?r 1 ?r), drop(?b ?g ?r)) ϕ 2 = (pick(?b ?g ?r 1 ), move(?r 1 ?r), drop(?b ?g ?r)) ϕ 3 = (goto(?r 1 ), pick(?b ?g ?r 1 ), move(?r 1 ?r), drop(?b ?g ?r)) ϕ 4 = (goto(?r 1 ), pick(?b ?g ?r 1 ), goto(?r), drop(?b ?g ?r))

. . . Some of these decomposition are redundant. To facilitate proof reading we want a more compact description of the HTN Methods. More precisely, we want minimizing the set of methods allowing to decompose observed compound tasks. Then, the HTN Methods learning problem can be reduce to a variant of the set cover problem [START_REF] Karp | Reducibility Among Combinatorial Problems[END_REF]) which is NP-Complete.

The Greedy algorithm [START_REF] Jungnickel | The Greedy Algorithm[END_REF]) is a classical way to approximate the solution in a polynomial time. At each stage of the Greedy algorithm, HierAMLSI chooses The main drawback of this approach is that it does not take into account the fact that there are dependencies between tasks. For instance, the optimal way to decompose the compound task move1ball(?b ?r) is ϕ 4 = (goto(?r 1 ), pick(?b ?g ?r 1 ), goto(?r), drop(?b ?g ?r)) (see Figure 3). So the compound task move1ball(?b ?r) depends of the compound task goto(?r). So, as long as all methods for the compound task goto(?r) have been generated, the Greedy algorithm will always prioritize ϕ 1 and ϕ 2 to ϕ 4 . Indeed, the decomposition ϕ 4 can be added to the current solution of the Greedy algorithm if and only if all methods for the compound task goto(?r) have been added.

Heuristic Approach

We propose a sound, complete and polynomial Heuristic approach (see Algorithm 1) taking into account dependencies between tasks. Figure 4 gives an example for the Gripper domain.

HierAMLSI starts by initializing the set of HTN methods using the decomposition ϕ observed during the observation generation step (see Section 4.1). For each Compound Task we have therefore a set of HTN Methods containing only Primitive Tasks and no Compound Task dependencies. For instance, for the Compound Task move1ball(?b ?r) we have the two following decomposition: ϕ 1 = (pick(?b ?g ?r 1 ), move(?r 1 ?r), drop(?b ?g ?r)) ϕ 2 = (move(?r ? ?r 1 )), pick(?b ?g ?r 1 ), move(?r 1 ?r), drop(?b ?g ?r))

Then, at each iteration AMLSI use the greedy algorithm to compute a new set of HTN Methods with an additional Compound Task dependency. Finally, if the new HTN Methods set is smaller than the one learned in the previous iteration, then it is retained. For instance, for the Gripper domain, suppose we have the two following decomposition for the Compound Task goto(?r):

ϕ 1 = () ϕ 2 = (move(?r 1 ?r)))
Then, the Greedy algorithm return only one decomposition for the Compound Task move1ball(?b?r): (goto(?r 1 )), pick(?b ?g ?r 1 ), goto(?r), drop(?b ?g ?r)).

Lemma 1. The Heuristic approach is sound and complete. The heuristic approach generates a set of HTN Methods M able to decompose all observed compound tasks in the dataset Ω.

Proof. During the observation generation step (see Section 4.1), for each generated Compound Task t, we have its final decomposition ϕ. So at the initialization step of the Heuristic approach, there at least one method able to decompose each observed Task. The initialization is therefore sound and complete. Moreover the following steps of the Heuristic approach generates methods decomposing as many tasks as the previous steps, then the Heuristic approach is sound and complete.

Lemma 2. The Heuristic approach is polynomial.

goto(?r) : ϕ 1 = (), ϕ 2 = move(?r 1 ?r) move2balls(?b 1 ?b 2 ?r) : ϕ 1 = (pick(?b 1 ?g 1 ?r 1 ), pick(?b 1 ?g 1 ?r 1 ), move(?r 1 ?r), drop(?b 1 ?g 1 ?r), drop(?b 1 ?g 1 ?r))

ϕ 2 = (move(?r ?r 1 ), pick(?b 1 ?g 1 ?r 1 ), pick(?b 1 ?g 1 ?r 1 ), move(?r 1 ?r), drop(?b 1 ?g 1 ?r), drop(?b 1 ?g 1 ?r)) move1ball(?b ?r) :

ϕ 1 = (pick(?b ?g ?r 1 ), move(?r 1 ?r), drop(?b ?g ?r)) ϕ 2 = (move(?r ? ?r 1 )), pick(?b ?g ?r 1 ), move(?r 1 ?r), drop(?b ?g ?r))

(a)

Step 0: Initialization with no compound task dependency goto(?r) : ϕ 1 = (), ϕ 2 = move(?r 1 ?r) move2balls(?b 1 ?b 2 ?r) : ϕ 1 = (goto(?r 1 ), pick(?b 1 ?g 1 ?r 1 ), pick(?b 1 ?g 1 ?r 1 ), goto(?r), drop(?b 1 ?g 1 ?r), drop(?b 1 ?g 1 ?r)) move1ball(?b ?r) :

ϕ 1 = (goto(?r 1 ), pick(?b ?g ?r 1 ), goto(?r), drop(?b ?g ?r))}

(b)

Step 1: Induction with 1 compound task dependence goto(?r) : ϕ 1 = (), ϕ 2 = move(?r 1 ?r) move2balls(?b 1 ?b 2 ?r) : ϕ 1 = (goto(?r 1 ), pick(?b 1 ?g 1 ?r 1 ), pick(?b 1 ?g 1 ?r 1 ), goto(?r), drop(?b 1 ?g 1 ?r), drop(?b 1 ?g 1 ?r)) move1ball(?b ?r) :

ϕ 1 = (goto(?r 1 ), pick(?b ?g ?r 1 ), goto(?r), drop(?b ?g ?r))}

(c)

Step n: Induction with n compound task dependence Table 1: Benchmark domain characteristics. From left to right, the number of Primitive Tasks, the number of Compound Tasks, the number of Methods and the number of Predicates for each IPC domain.

Proof. First of all, we have O(|I + |) 2 nodes in the DFA. Then, in the worst case, we have a possible HTN Method for each node pair, then we have O(|I + | 2 ) possible HTN Methods in the DFA. Then, the complexity of the Greedy algorithm is O(|I + | 3 ) in term of tested decomposition. Finally, according to the algorithm 1, the Greedy algorithm is repeated |C| 2 times. Finally, the complexity of the Heuristic approach is O(|C| 2 .|I + | 3 ).

Exprimentation

The purpose of this evaluation is to evaluate the performance of HierAMLSI though two variants: (1) we evaluate the performance of HierAMLSI when only

Experimental Setup

Our experiments are based on 5 HDDL [START_REF] Höller | HDDL: An Extension to PDDL for Expressing Hierarchical Planning Problems[END_REF][START_REF] Höller | Hierarchical Planning in the IPC[END_REF]) domains (see Table 1) from the IPC 2020 competition: Blocksworld, Childsnack, Transport, Zenotravel and Gripper. HierAMLSI learns HTN domains from one instance. To avoid performances being biased by the initial state, Hier-AMLSI is evaluated with different instances. Also, for each instance, to avoid performances being biased by the generated observations, experiments are repeated five times. All tests were performed on an Ubuntu 14.04 server with a multi-core Intel Xeon CPU E5-2630 clocked at 2.30 GHz with 16GB of memory. PDDL4J library [START_REF] Pellier | PDDL4J: a planning domain description library for java[END_REF] was used to generate the benchmark data.

Evaluation Metrics

HierAMLSI is evaluated using the accuracy (Zhuo, Nguyen, and Kambhampati 2013) that measures the learned domain performance to solve new problems.

Formally, the accuracy Acc = N N * is the ratio between N , the number of correctly solved problems with the learned domain, and N * , the total number of problems to solve. In the rest of this section the accuracy is computed over 20 problems. The problems are solved with the TFD planner [START_REF] Pellier | Totally and Partially Ordered Hierarchical Planners in PDDL4J Library[END_REF] provided by the PDDL4J library. Plan validation is done with VAL, the IPC competition validation tool [START_REF] Howey | VAL's progress: The automatic validation tool for PDDL2. 1 used in the international planning competition[END_REF].

Discussion

Table 5 shows the accuracy of HierAMLSI obtained on the 5 domains of our benchmarks when varying the training data First of all, we observe that when HierAMLSI learns only the set of HTN Methods, learned domains are generally optimal (Accuracy = 100%) with 600 tasks whatever the experimental scenario. Also, 100 tasks are generally sufficient to learn accurate domains (Accuracy > 50%). Then, when HierAMLSI learns both action model and HTN Methods performances are similar when observations are noiseless. However, when observations are noisy, performances are downgraded for some domains: Blocksworld and Childsnack when observations are complete and Blocksworld, Transport and Childsnack when observations are partial. This is due to the fact that there are learning error in the action model learned. However, learned domains remain accurate when there are at least 300 tasks in the training dataset.

To conclude, we have shown experimentally that Hier-AMLSI learns accurate domains. More precisely, when the action model is known, HierAMLSI generally learns optimal domains. Also the performances are downgraded when HierAMLSI has to learn the action model in addition to the set of HTN methods, but the learned domains remain accurate. Performance degradation are due to the learning errors in the action model.

Related Works

Many approaches have been proposed to learn HDDL domains. These approaches can be classified according to the input data of the learning process. The input data can be plan "traces" obtained by resolving a set of planning problems, annoted plans (see Figure 6a), decomposition tree (see Figure 6b) or random walks. The input data can contain in addition to the tasks also states which can be fully observable (FO) or partially observable (PO), or noisy. Also, these approaches can be classified according to the output. The output can be the action model, the set of HTN Methods and HTN Methods preconditions. Table 2 summarises these classifications.

A first group of approaches only learns the set of HTN Methods. First of all, [START_REF] Xiaoa | Learning HTN Methods with Preference from HTN Planning Instances[END_REF]) takes as input a set of plan traces and HTN Methods and proposes an algorithm to update incomplete HTN Methods by task insertions. Then HTN-Maker [START_REF] Hogg | HTN-MAKER: Learning HTNs with Minimal Additional Knowledge Engineering Required[END_REF] and HTN-Maker ND [START_REF] Hogg | Learning hierarchical task networks for nondeterministic planning domains[END_REF] takes as input plan trace generated from STRIPS planner and annoted task provided by a domain expert. Then, [START_REF] Hogg | Learning methods to generate good plans: Integrating htn learning and reinforcement learning[END_REF] proposed an algorithm based on reinforcement learning. Then, [START_REF] Li | Learning probabilistic hierarchical task networks as probabilistic context-free grammars to capture user preferences[END_REF]) proposed an algorithm taking as input only plan traces. This algorithm builds, from plan traces, a context free grammar (CFG) allowing to regenerate all plans. Then, methods are generated using CFG: one method for each production rule in the CFG. Then [START_REF] Garland | Learning hierarchical task models by demonstration[END_REF] and [START_REF] Nargesian | LHTNDT: Learn HTN Method Preconditions using Decision Tree[END_REF] proposed to learn HTN Methods from annoted plan. Annoted plan are plan segmented with the different tasks solved. Figure 6a gives an example of annoted plan. However, obtaining these annotated examples is difficult and needs a lot of human effort.

A second group of approach learns HTN Methods preconditions. First of all, the CAMEL algorithm [START_REF] Ilghami | CaMeL: Learning Method Preconditions for HTN Planning[END_REF] learns HTN Methods and their the preconditions of HTN Methods from observations of plan traces, using the version space algorithm. An annoted task is an triplet (n, P re, Ef f ) where n is a task, P re is a set of proposition known as the preconditions and Ef f is a set of atoms known as the effects. These approach use annoted task to build incrementally HTN Methods with preconditions. Only HTN-Learner proposes to learn Action Model and HTN Methods from decomposition trees. A decomposition tree is a tree corresponding to the decomposition of a method. Figure 6b gives an example of decomposition tree.

Conclusion

In this paper we have presented HierAMLSI, a novel algorithm to learn HDDL domains. HierAMLSI is built on the AMLSI approach. HierAMLSI is composed of four steps. The first step consists, as AMLSI, in building two training sets of feasible and infeasible action sequences. In the second step, HierAMLSI induces a DFA. The third step is the generation of the HTN Methods, and the last step learns HTN Methods preconditions and the action model. Our experimental results show that HierAMLSI is able to learn accurately both action models and HTN Methods from partial and noisy datasets.

Future works will focus on extending HierAMLSI in order to learn more expressive action model.

  move2balls(b 1 b 2 r b ) has been decomposed by primitive tasks (move(r a r b ), pick(b 1 r r b ), pick(b 2 l r b ), move(r b r a ), drop(b 1 r r a ), drop(b 2 l r a )) in node 0 and reached the node 6. Then we add the following transitions in the DFA γ(move2balls(b 1 b 2 r b ), 0) → 6. Figure 2 gives an example of Task DFA.

Figure 2 :Figure 3 :

 23 Figure 2: DFA learning steps The Primitive Task DFA is the DFA containing only Primitive Task, i.e. black transition, and the Task DFA contains Compound Tasks, i.e. red transitions, in addition

Figure 4 :

 4 Figure 4: HTN Methods learning example

Figure 5 :

 5 Figure 5: Performance of HierAMLSI when the training data set size increases in number of tasks in terms of Accuracy.

Table 2 :

 2 State-of-the-art action model learning algorithms. From left to right: the kind of input data, the kind of environment: Fully Observable or Partially Observable, the maximum level of noise in observations, the kind of output Then, the HDL algorithm[START_REF] Ilghami | Learning to Do HTN Planning[END_REF] takes as input plan traces. For each decomposition in plan traces, HDL checks if there exist a method responsible of this decomposition. If not, HDL creates a new method and initializes a new version space to capture its preconditions. Preconditions are learned in the same way as in the CAMEL algorithm.

	Algorithm	Input	Input Environment Noise Level Action Model HTN Methods HTN Methods Preconditions Output
	(Xiaoa et al. 2019)	Plans	FO	0%		•
	HTN-Maker	Plans	FO	0%		•
	(Hogg, Kuter, and Munoz-Avila 2010)	Plans	FO	0%		•
	HTN-Maker ND	Plans	FO	0%		•
	(Li et al. 2014)	Plans	FO	0%		•
	(Garland and Lesh 2003)	Annoted Plans	PO	0%		•
	LHTNDT	Annoted Plans	FO	0%		•
	CAMEL	Plans	FO	0%		•	•
	HDL	Plans	FO	0%		•	•
	HTN-Learner	Decomposition Trees	PO	0%	•	•	•
	HierAMLSI	Random Walks	PO	20%	•	•	•

Acknowledgments

This research is supported by the French National Research Agency under the "Investissements d'avenir" program (ANR-15-IDEX-02) on behalf of the Cross Disciplinary Program CIRCULAR.