D Van Zwynsvoorde

T Simeon

R Alami

Incremental Topological Modeling using Local Vorono -like Graphs

In the eld of mobile robotics, one important issue is to allow the robot to navigate in an a priori unknown and non speci c large scale environment. Large dimensions raise strong limitations of geometric modeling, and topological or mixed metric-topological models are now studied to better t the problem.

In this paper, we present a new method for incrementally building a topological model of an indoor environment from sensor range data. The approach consists in merging each local perception of the topology with the current state of the global graph. This local topology is captured through the construction of a Vorono -like graph that takes into account not only visible features but also visibility constraints (hidden regions, limited sensing ranges, ...). We give the outline of the method and show rst encouraging results on real data.

Introduction

This paper addresses the issue of long term navigation in an a priori unknown, large scale environment. Large dimensions have dramatic consequences on localization since odometric pose estimate increasingly degrades as the robot moves around.

There exist some e ective simultaneous mapbuilding and localization methods (e.g 2, 3]) but they su er from a rather costly learning scheme in terms of memory space and computation time. Indeed, they are based on localizing the robot in a global metric map of the environment, which becomes di cult in a large scale environment because of large position uncertainty. Some authors also de ned sensor based strategies to deal with uncertainties 1, 4, 5, 6] but they assume some a priori knowledge and/or are limited to local navigation.

In such wide environments, a topological representation 8, 10, 11, 13, 14] may be more e ective because it provides a compact model that allows fast route planning and reduces the amount of stored data, as well as the computation time for learning. There are two common approaches to topological modeling. Thrun 13] proposes a strategy that constructs a topological graph over a metric representation of the environment by partitioning the map into adjacent regions separated by narrow passages. It gains from both types of representation (metric and topological) but also su ers from the disadvantages of the metric map. Navigation techniques using such topological/metric descriptions have also been investigated (eg. 14] based onto Markov processes).

The other class of methods tries to construct a topological model directly from sensor data. [START_REF] Kuipers | A robot Exploration and Mapping Strategy Based on a Semantic Hierarchy of Spatial Representations[END_REF], Dubek proposes a hybrid topological-metric model where nodes of the topological graph are local metric descriptions of \Islands of reliability". (see also 12] for a similar approach). Choset describes in 8] a sensor-based exploration procedure allowing a mobile robot to trace the Vorono of the environment from sonar inputs. The system incrementally constructs a graph connecting meet points lying at the intersection of Vorono 's edges. At each step, sonar data are used to maintain the robot onto the graph and to locally extend the edge currently traced with a single point.

Our topological model also derives from the Vorono structure. However, the approach we propose to build this model exploits the richness of the information that can be collected by a laser range nder. Local Voronolike graphs are computed from the segmented range data and used to incrementally build the topological model. The potential interest of such approach is that matching local graphs with the partially constructed graph may allow to better maintain the reliability of the topological map.

The paper is organized as follows. Section 2 gives an overview of the approach. Section 3 presents how the local model is obtained from sensor data. Section 4 proposes a procedure for local merging of local and global models, along with a loop detection scheme to allow modeling cyclic environments. Section 5 gives some preliminary results on a real environment.

2 Overview of the Approach

The problem addressed in this paper is illustrated by gure 1. The left part shows the whole set of segments constructed from the sensor range data collected by a mobile robot during an experiment where the robot moved within the lab. One can note in the gure that many segments are duplicated and seem as shifted and twisted. This is the consequence of the dead-reckoning error accumulated by the robot during its displacements. The aim of this work is to be able of incrementally constructing a topological model of the map, such as the one depicted at the right part of the gure. The approach is based on the assumption that the explored world can be reliably modeled by segments obtained from a range nder performing an horizontal scan. We consider a laser sensor that has both limited distance and angular ranges (see Fig 3). Given a set of points collected by the sensor, a segmentation algorithm produces a sequence of possibly disconnected segments modeling the sensed obstacles. Discontinuities in the sequences of segments may be due to occlusions, specular limitation of the sensor (see Fig 4), or \unperfect" segmentation. By appending escapelines departing from each discontinuity, we construct a generalized star-shaped visibility polygon which represents the current free space area perceived by the robot.

Line Segments

Global Topological Graph

Occlusion area

Local Vorono Diagram

The Vorono Diagram is a classical structure for path planing 1] de ned by the set of points equidistant to at least two obstacles. Three-equidistant points correspond to possible bifurcations. Thus the Vorono Diagram provides a natural way to de ne a graph that explicitly captures the topology of the free space.

Considering the generalized visibility polygon as the free space area, we extract the local topology using an algorithm originally presented in 7]. The principle consists in tracing incrementally the Vorono by moving a virtual point denoted R, initially placed at the sensor's position. An initialization step rst computes the retraction of R onto the diagram. Successive steps computed based on the geometric analysis of the neighborhood of R (vertices on the boundary, nearest obstacles...) trace the current edge until R reaches a 3-equidistant point. A node is created at this position, with departing edges. The algorithm then selects one unexplored edge of the new node and continues tracing. It terminates when all the edges have been explored. As illustrated by Fig 5, each of the portions of the edges is labeled as de nitive (if based only on boundary segments) or temporary ( when at least one escapeline is involved). When one temporary portion intersects with another edge, the corresponding node is also labeled as temporary. Temporary nodes are likely to drift around their current position inside a presencearea. This area can be over-estimated by a circular region de ned by the extremal positions of the implicated escape-lines (see Fig. 6). This area will eventually contain at least one de nitive node.

Structure of the Local Map

Mainly because of unperfect (noisy) data segmentation, this rough graph must be ltered to keep the useful edges for navigation. In particular we discard the edges based on too at concave vertices or due to narrow passages that do not allow the robot to navigate.

Nodes of the graph are described by their uncertain position (because of the uncertain position of the robot), the anticlockwise ordered set of edges emanating from them (see Fig 7) indexed with relative departing angles, a temporary/de nitive attribute, and a con dence degree that re ects the persistence of the node along the successive local views (see next section).

Each of the edges is described by a sequence of portions. A portion is either a straight or curve segment which has an approximate length. It is also assigned a temporary/de nitive attribute and distance to the environment is stored for each extremity. Global attributes of the entire edge can be derived from a simple integration of the local attributes of the successive portions. At each step of the exploration, the system gets the topology of the nearby environment. We assume that we can get a retraction of the robot position onto the graph and determine which is the direction of displace-ments on the corresponding edge. Consequently, we naturally use the local graph to locally update the ongoing part of the current global graph.

Moreover, considering small displacements, we expect that an important part of the local graph will remain qualitatively identical over successive views. However, temporary nodes will have slightly changed and some nodes may have appeared or disappeared, as illustrated in gure 8. The gure shows the local graph constructed at two close positions of the robot. In the right gure, the specular phenomenon present in the left one has disappeared so that two new nodes can be modeled, while one node already present in the left gure has turned from temporary to de nitive.

We next describe a procedure allowing to identify these local changes. Then we address the issue of loop detection to handle situations where the robot discovers again a previously modeled region of the environment.

Presence Area

Local Growth of the Global graph

The point is to determine, for each node N l of the local graph if it is already modeled in the global graph or not. The key idea is that the topological feature represented by N l in the local graph is already modeled in the global one if we can reach a node N g in the global graph by mean of a path topologically equivalent to the one that leads to N l in the local graph. This accessibility procedure works as follows. Note D the closer extremity of the local graph to the robot position. We construct a path from D to N l as a set of nodes fL i g to be traversed. Now examine this path \node-by-node". To switch from one node L i to the next, we de ne manoeuvres. Arriving at Continuing this accessibility procedure for each L i leads to one of the following situations :

Both N l and a corresponding global node N g are reached. Then N l is used to update N g 's description, and N g 's con dence degree is increased (eg. nodes L 2 and L 3 of g 10). N l is reached but a still pending edge of the global graph has been reached. N l is connected to this edge with a default value of the con dence degree (eg. node L 4 ). N l is situated before the presence area of the corresponding node reached in the global graph. This can occur when the perceptual situation has changed so that a new node is discovered. N l is inserted in the global graph with a default condence degree (eg. node L 1 ). One L i lies after the presence area of the reached global node G i . Then we determine if it is consistent with the new robot position, in particular considering possible occlusions and specularity. If it is, the accessibility procedure continues without updating anything; if not, the con dence degree of G i is decreased and a copy of G i is inserted in the local graph to allow continuing the procedure. N l stands inside the presence area of the reached global node G i , but they have di erent topologies. This is possible only if G i is constructed over escape-lines that have strongly evolved in the local perception (see g 8). N l is inserted but with a low degree of con dence if it is not the case. This accessibility procedure performs a local growth of the global model, without considering closing-loop situations.

Loop Detection

At each step of the construction, the system must detect if the robot comes to explore again a region it has already modeled or not. Loop detection is a tricky issue of topological modeling. In large scale environments, large dead-reckoning errors prevent from identifying the nodes using their position. To perform node matching, some authors have proposed various criteria based on the nodes description. In 8] Choset uses four stable characteristics of the nodes. First, distance to nearest obstacles is examined. The number of emanating edges is also a relevant feature. A more distinctive criterion is the comparison on departure angles of the emanating edges. The last useful characteristic of a node is the distinction between edges terminating on obstacles boundary (named blind edges) and others.

In 9], the geometrical description of the edges is also used, but as the system builds the edges \pointby-point", it must trace the edge to get su cient information. In our model we already have the description of the edges (or at least part of it) when trying to identify a node, so we can use comparison of edges as a more discriminating criterion, in addition to the previous criteria.

The matching of edge description consists in comparing edge curvature and distance to obstacles all along the common de nitive part of both edges. Note that all necessary geometrical informations are stored into the edge structure of the graph.

For each de nitive node in the local graph, applying this multi-criteria matching may result in several candidate nodes if the robot position (and consequently nodes position ) is too uncertain. Then we use a classical scheme for generation and validation of hypothe- Let us explain this procedure on a (not so) simple example. Figure 11 shows a situation where the robot discovers a node (node 8) topologically similar to several candidate nodes (nodes 1, 2 and 3) already present in the model. Since the robot has covered a large distance, its position uncertainty has become important enough so that nodes 1, 2 and 3 lie in the uncertain domain of node 8. Then a list of hypotheses is initiated. Figure 11 shows how new nodes allow validation of some of the hypotheses. Node 9 is supposed to be one of nodes 1,2 or 3. Two hypotheses are made that respectively connect node 1 to node 2 and node 2 to node 3. The procedure continues until node 11 leads to the conclusion that none of the hypotheses can be validated because node 5 is topologically di erent from node 4 which is the only following node from node 3. Thus nodes 8, 9 and 10 are new ones. Note that even if there is only one candidate node at the rst step of a closing loop situation, the loop is not necessary e ective. So the hypothesis must be stable over a few steps to be con rmed (or not). One may also note that the local topogical map computed with our approach may allow a more e ective matching since the topology of the neighbor nodes can also be used when matching local nodes into the gobal model. Such neighboring information possibly allows to eliminate some potential candidates , leading to a reduced number of hypotheses. This is illustrated by the slightly modi ed example of Fig. 12. Here, the robot perceives two nodes (nodes 8 and 9) when arriving in the vicinity of node 8, and thus can directly conclude that node 8 cannot be either of nodes 1, 2 nor 3 (because none of these nodes is connected to a node such as node 9).

When a loop has been identi ed, we make a systematic comparison between each encountered node respectively in the old and new part of the hypothesis, backtracking from the last node of the hypothesis until the end of the revisited region is detected. Then the loop is connected erasing the revisited part of the model, and the robot position is corrected.

Experimental Results

We performed some experiments using our robot Hi-lare2 as a mobile sensor to collect laser range data. The robot was manually driven across the lab and the segmented data were stored every 20 cm. The gures below detail some steps of the construction of the topological graph shown in Fig. 1. Each of the local view contains an of 20 segments in a room-like neighborhood, and less than 10 segments in a corridor. In the rst case the average computation time of the local graph is about 0.5s. In the corridor case, the local graph is obtained after less than 20ms.

One can also note over the successive steps that the global graph is roughly divided into a de nitive and a temporary part. Recently discovered nodes are often temporary but they are updated in further exploration so that the temporary part of the global graph tends to converge towards the real skeleton of the environment. Note in gure (b) that because of occlusions and specularity, the local graph only builds a \very temporary" node (at the bottom of the gure), whereas gure (c) shows that the following steps allowed modeling more nodes. Figure (d) explicitly shows that the local graph may not be connex. The last gure (e) corresponds to the robot entering a corridor. Note that despite a quite large sensing range (over 10 meters) the local model (side walls) remains rather limited because of specularity.

Future Work and Conclusion

We have developed a new approach for directly constructing a topological model of an indoor structured environment from local range data. The model is obtained incrementally by merging each local topological graph with the current state of the model. Some encouraging results have been obtained on a real environment and we are currently validating the loop detection scheme.

Future work will concentrate on extracting exploration and navigation strategies, and also on demonstrating the robustness of the approach for long range navigation.

Figure 1 :

 1 Figure 1: Topological map built in a real environment

FigureFigure 2 Figure 3 :

 23 Figure 2: ApproachFigure 2 summarizes the scheme we adopted for topological modeling. As the robot is exploring the

Figure 4 :

 4 Figure 4: Escape-lines: Thick solid lines are the segments built from sensor data and thin solid lines are the escape-lines.

Figure 5 :

 5 Figure 5: Local Vorono graph. The graph drawn with thick solid lines represents the de nitive part of the local model, while dashed edges are based on escape-lines and are only temporary. Empty squares are placed over each node.

Figure 6 :

 6 Figure 6: Presence area for one node: N' and N" are the equivalent position of the node N without considering each of the two escape-lines.

  Figure 7: Topological Node

Figure 8 :

 8 Figure 8: Local evolution: lled black squares represent de nitive nodes whereas empty squares denote temporary nodes. From the left to the right gure the specularity disappeared and a new node is created lying inside the presence area of an existing node

Figure 9 :

 9 Figure 9: Manoeuvres arriving at a node one node, the robot can either keep forward (KF), turn right (TR), turn left (TL), avoid right (AR) or avoid left (AL) (see gure 9). Note D 0 the retraction of D onto the global graph. Starting from D 0 we examine the rst encountered node G 0 and compare it with the rst local node L 0 . For the equivalence to be acceptable, both nodes L 0 and G 0 must o er compatible manoeuvres.

Figure 10 :

 10 Figure 10: Local evolution: Thick lines are the global model and thin lines are the local graph.

Figure 11 :

 11 Figure 11: Loop situation and generated hypotheses ses.Let us explain this procedure on a (not so) simple example. Figure11shows a situation where the robot discovers a node (node 8) topologically similar to several candidate nodes (nodes 1, 2 and 3) already present in the model. Since the robot has covered a large distance, its position uncertainty has become important enough so that nodes 1, 2 and 3 lie in the uncertain domain of node 8. Then a list of hypotheses is initiated. Figure11shows how new nodes allow validation of some of the hypotheses. Node 9 is supposed to be one of nodes 1,2 or 3. Two hypotheses are made that respectively connect node 1 to node 2 and node 2 to node 3. The procedure continues until node 11 leads to the conclusion that none of the hypotheses can be validated because node 5 is topologically di erent from node 4 which is the only following node from node 3. Thus nodes 8, 9 and 10 are new ones. Note that even if there is only one candidate node at the rst step of a closing loop situation, the loop is not necessary e ective. So the hypothesis must be stable over a few steps to be con rmed (or not).

Figure 12 :

 12 Figure 12: Loop detection using neighboring information

Figure 13 :

 13 Figure 13: Experiment: The whole set of sensed segments and the rst ltered local graph The dimensions of the environment explored during this experiment are about 20 10 meters and the topological graph of Fig. 1 contains 22 nodes (10 of them are de nitive nodes). The initial state of the global graph (corresponding to the rst local graph) is shown in Fig. 13. This local graph was obtained after ltering of the graph displayed in Fig. 14.

Figure 14 :

 14 Figure 14: First local graph before ltering

  A few steps of the global constructionFor each step, the left part corresponds to the local graph displayed with the local data (thick solid lines are the boundary segments and thin dashed lines are the escape-lines). The right part of the gure shows the new global graph after graph matching. Black squares denote de nitive nodes whereas grey ones are only temporary. Thick solid lines are de nitive edges of the graph and thin dashed ones represent temporary portions.

Table ,

 ,

	Chairs & Cupboard	Office entry


Vibrotactile thresholds on the wrist for vibrations normal to the skin Quentin Consigny, Nathan Ouvrai, Arthur Paté, Claudia Fritz and Jean-Loïc Le Carrou

Abstract-In order to facilitate the development of wrist-worn vibrotactile devices, detailed knowledge about how vibrations are perceived by the users is needed. In particular, perceptual thresholds in amplitude are really important. Thresholds have been measured in the literature for other areas of the body, but given the variability reported between areas (shape of the threshold curve, position of maximum sensitivity), thresholds on the wrist can not be inferred from previous measurements and must be measured. The amplitude thresholds for vibrations normal to the skin surface were evaluated on 28 participants, with a three interval forced choice method. They were measured for 7 frequencies that are classical in the literature about vibrotactile perception [START_REF] Salzer | Vibrotactile "onthigh" alerting system in the cockpit[END_REF]40,80,160,250,320,and 640 Hz). The classical U-shape of the amplitude-threshold curve is observed, with a maximum sensitivity at around 160 Hz, which differs from other body areas, but confirms recent results obtained for vibrations parallel to the skin surface of the same body area. The sensitivity thresholds of vibrotactile signals appear to be in the micrometer range.

I. INTRODUCTION

The number of vibrotactile devices worn on the wrist has increased in recent years. These devices can be designed for research activities [START_REF] Ma | Vibrotactile wristband for warning and guiding in automated vehicles[END_REF]- [START_REF] Paredes | CalmMeNow[END_REF] or developed for a "general public" use (connected watch or haptic metronome 1 for example). Wrist-mounted vibrotactile devices have three advantages:

(1) a good sensitivity due to a relatively high density of mechanoreceptors (there are 10 times more receptors on the wrist than on the back) [START_REF] Corniani | Tactile innervation densities across the whole body[END_REF]; (2) a low bone conduction of sound [START_REF] Lenhardt | Bone-conduction propagation in the human body: implications for high-frequency therapy[END_REF] so that purely vibrotactile applications can be targeted, where the vibrotactile message can be broadcasted in noisy environments and accessible to hearing-impaired users and (3) a good accessibility, allowing to keep the hands free. The emergence of new vibrators (electro-dynamic transducers dedicated to vibrotactile applications) allows diversifying the offer and the capacities of these bracelets [START_REF] Choi | Vibrotactile display: Perception, technology, and applications[END_REF]. The miniaturization and democratization of the technologies used in the generation of vibrotactile signals (sound card, transducer, battery, etc.) facilitates the design of wrist-worn vibrotactile devices. In recent years, new technology industries have developed new applications based on the use of the sense of touch. Users report a real interest in this kind of technology. Vibrotactile wristbands could therefore continue to develop in the coming years.

This work was funded by a grant from the French National Research Agency (ANR) as part of the "Staccato" Project (ANR-19-CE38-0008-01). This project has been validated by the research ethics committee of the University of Paris. IRB NUMBER : 00012020-47. The research team thanks the 28 people who participated in this experiment and Hugues Genevois for the valuable discussions. 1 Pulse made by Soundbrenner.

For this, it is necessary to have a good knowledge of the users' feelings in order to design the vibrators and to anticipate the perception of the proposed vibrotactile signals. A first approach, classically used in psychophysics, is to analyze, both in amplitude and frequency, the user's sensitivity to vibrotactile signals. Detection of these signals is an essential parameter that is important to characterize, in particular through the identification of perceptual thresholds [START_REF] Verrillo | Investigation of some parameters of the cutaneous threshold for vibration[END_REF].

Nevertheless, the perception and detection of vibrotactile signals is dependent on the context of the experiment, such as the type of vibrator [START_REF] Xie | A review of smart materials in tactile actuators for information delivery[END_REF], [START_REF] Vorobyov | Comparison of 3d printed vibro-tactile actuator with permanent magnet only and standard vibro-actuator for prosthetic feedback devices[END_REF], the force exerted on it [START_REF] Papetti | Vibrotactile sensitivity in active finger pressing[END_REF]- [START_REF] Pra | Effects of vibration direction and pressing force on finger vibrotactile perception and force control[END_REF], the shape in contact with the skin [START_REF] Verrillo | Effect of contactor area on the vibrotactile threshold[END_REF], [START_REF] Schmidt | Larger contactor area increases low-frequency vibratory sensitivity in hairy skin[END_REF] and the user themselves (hairiness [START_REF] Schmidt | Larger contactor area increases low-frequency vibratory sensitivity in hairy skin[END_REF], [START_REF] Verrillo | Vibrotactile thresholds for hairy skin[END_REF], age [START_REF] Seah | Normal values for thermotactile and vibrotactile thresholds in males and females[END_REF], temperature [START_REF] Harazin | Effect of changes in finger skin temperature on vibrotactile perception threshold[END_REF]). Naturally, the area of the body that is excited induces significant variations: this was measured on the finger [START_REF] Verrillo | Vibration sensation in humans[END_REF]- [START_REF] Birnbaum | A systematic approach to musical vibrotactile feedback[END_REF], the arm [START_REF] Mirzaei | Effects of using vibrotactile feedback on sound localization by deaf and hard-of-hearing people in virtual environments[END_REF], the back [START_REF] Mirzaei | Effects of using vibrotactile feedback on sound localization by deaf and hard-of-hearing people in virtual environments[END_REF], the thigh [START_REF] Salzer | Vibrotactile "onthigh" alerting system in the cockpit[END_REF], the skull [START_REF] Kim | Rendering vibrotactile flow on backside of the head: Initial study[END_REF] or the mouth [START_REF] Barlow | Adaptive vibrotactile threshold estimation of the glabrous hand and perioral face following MCA stroke[END_REF]. For the wrist, detection thresholds were determined by exciting the skin parallel to its surface (the motion of the vibrator is along the skin, mostly using shear forces) [START_REF] Aevarsson | Vibrotactile threshold measurements at the wrist using parallel vibration actuators[END_REF]. However, some vibrotactile bracelets (especially those using voice coils), deform the skin in the direction normal to this surface. Therefore, thresholds measurements with vibrations normal to the skin using such devices are needed and would complete the work previously done.

In this article, amplitude thresholds for vibrotactile signals were measured with an electrodynamic transducer, placed on the internal face of the wrist. In parallel, a vibratory characterization of the vibrators was carried out, with the aim of identifying the frequency response, sensitivity and distortion behavior of the vibrators. After a description of the method and the experimental protocol (Section II), the thresholds of vibrotactile perception are presented (Section III) and discussed (Section IV) before concluding (Section V).

II. METHOD

A. Material

The vibrotactile device used in this work is made of an electrodynamic vibrator2 maintained in contact on the volar wrist3 by a watchband 4 which can be adjusted. This model is powerful enough to generate easily perceptible signals while not exceeding the size of a watch (a size that is accepted by a large part of the population). The fixed part of this vibrator is inserted in a wooden plate of 5 mm thickness. The mobile part is made of a crown of a few mm 2 on which a Plexiglas disk (3 cm diameter, 3 mm thickness) is glued to facilitate the perception of vibrations [START_REF] Verrillo | Psychophysics of vibrotactile stimulation[END_REF], [START_REF] Verrillo | Sensation magnitude of vibrotactile stimuli[END_REF]. Studies on vibrotactile thresholds on the finger show that a decrease of the area of skin in contact with the vibrator lowers the perceptual thresholds and their evolution with frequency. For areas larger than a few tenths of a cm 2 , the effect dissipates. This plate (≈ 7 cm 2 ) thus makes it possible to get rid of this effect. The shape of the contactor (the Plexiglas disk) influences the perception threshold values [START_REF] Verrillo | Effect of contactor area on the vibrotactile threshold[END_REF]. The thickness of the disc prevents it from deforming under pressure. Thus, the shape of the contactor area is fixed, and is flat, neither concave, nor convex (see Figure 1). The notches of the watchband are a succession of 5mm wide magnets. The participants had to fit the bracelet by themselves so as the tightening felt comparable to the one used for a watch. Generally, once a participant found a comfortable tightness, adding a notch to the bracelet would make the vibrator lift off the skin while removing a notch generated too much tension on the bracelet so the magnets did not resist. With this tightening, participants exerted forces between 0.6 and 1.4 N (the mean value is 1 N, see Section VI), which is range of value already observed in a comparable study [START_REF] Pra | Effects of vibration direction and pressing force on finger vibrotactile perception and force control[END_REF].

The bracelet had to be worn on the wrist of the nonpreferred hand (i.e., left wrist for right-handed participant)

During the experiment, the participants were asked to keep the hand open and relaxed with the palm facing up. The participants' wrist rested on a piece of foam (see Figure 1). The participants wore earplugs 5 and ear muffs 6 in order not to hear the vibrator's noise (no participant reported having heard the stimuli). The earplugs and ear muffs each have an attenuation of 37 dB.

The vibrotactile signals were generated on a computer, then transmitted to a sound card 7 , then amplified by an amplifier 8 to finally be sent to the vibrator. During the characterization and the threshold experiment, the output volume of the computer, the gains of the sound card and the amplifier were fixed. The maximum voltage at the terminals of the transducer was 3.3 V (peak to peak). This setup is visible in Figure 1.

B. Experimental characterization of the vibrator

In order to perform perceptual threshold measurements, it is important to characterize the vibrator while being used on the skin, in terms of (1) its harmonic distortion (check whether additional/unwanted frequencies are generated by the system: if not, it ensures that only the wanted frequency is tested) and (2) its voltage-displacement law in order to compensate for the frequency response of the vibrator.

For the first characterization, the device is fixed on the volar wrist of a subject (see Figure 2). This measurement is long and tedious, so it was only done with one participant. For each frequency, bursts of 25 periods of a sine wave are sent to the vibrator. The amplitude of the first burst corresponds to the amplitude of the smallest threshold identified in Section III (see table II), decreased by 15 dB. Each burst has an amplitude 1 dB higher than the previous one, until the maximum possible amplitude (of 3.3 Volts peak to peak) is reached. The velocity of the vibrator is measured with a laser vibrometer 9 and its controller10 (see Figure 2). A Fourier transform is calculated on each burst. The harmonic distortion of the device is evaluated by computing the ratio between the sum of the amplitudes of the harmonics and the amplitude of the fundamental frequency [START_REF] Mayer | Investigation of vibrotactile transducers for a bone conduction sensory feedback system[END_REF]. This ratio is named THD (for Total Harmonic Distortion) and it is obtained as follows:

T HD = 100 × ∞ h=2 v 2 h v 2 1 , (1)
where v h is the RMS amplitude of the h-th harmonic. The THD is systematically measured for every generated burst.

Table II-B shows the highest THD values measured. These values are low enough to consider that, for a sinusoidal voltage input, the movement of the transducer is sinusoidal with the same frequency [START_REF] Mayer | Investigation of vibrotactile transducers for a bone conduction sensory feedback system[END_REF]. The same data is used for the second characterization, as the velocity can be integrated to estimate the instantaneous position of the vibrator and then its RMS value. The RMS position calculated on a period can be considered constant over time (we consider that a low THD indicates a low distortion of the waveform, and that this distortion evolves little in time). In order to limit the integration errors, these instantaneous and RMS positions are calculated on each period, the measurement is reproduced for the 15 periods of the burst and an average is calculated over these 15 measurements. Thus, the voltage-displacement law can be estimated (see Figure 3). The thresholds can therefore be expressed in displacement in order to facilitate the comparison with other published results.

C. Experimental protocol for the perceptual threshold

The amplitude threshold is estimated for 7 frequencies classically used in the literature [START_REF] Salzer | Vibrotactile "onthigh" alerting system in the cockpit[END_REF]40,80,160,250, 320 and 640 Hz) [START_REF] Verrillo | Effect of contactor area on the vibrotactile threshold[END_REF], [START_REF] Gescheider | A four-channel analysis of the tactile sensitivity of the fingertip: frequency selectivity, spatial summation, and temporal summation[END_REF] and covering most of the frequency range to which humans are sensitive via the vibrotactile modality [START_REF] Verrillo | Investigation of some parameters of the cutaneous threshold for vibration[END_REF], [START_REF] Verrillo | Sensation magnitude of vibrotactile stimuli[END_REF]. An "up-down" method associated with a "three interval forced choice" (3-IFC) procedure is used to determine the perceptual thresholds.

Figure 5 shows the threshold search method. For each frequency, a series of triads was presented to the participant (the order of these frequencies is randomized and balanced across participants). Each triad was composed of 3 one-second sequences (called intervals). Only one of those, chosen randomly, contained a stimulus. Between 2 intervals was a pause of one second. The participant could visualize on a screen (positioned in front of them) the current interval (see Figure 1). The workflow of the perceptual experiment is illustrated in Figure 4. At the end of the triad, the participant was asked to indicate in which interval they thought they perceived a stimulus. As long as the interval was correctly identified, a new triad at the same frequency was presented with a signal amplitude decreased by 3 dB. If not, the amplitude was increased by 3 dB and the process was repeated until the signal was identified 3 times consecutively. In a second step, the amplitudes were increased in 1 dB steps and 3 identifications in a row were systematically requested to consider the stimulus as perceived. This up-down pattern was repeated 3 times. The threshold is taken as the last amplitude. In order to reduce the duration of the experiment, the starting amplitude was set based on the results of a preliminary test (it was chosen between 3 and 26 dB below the maximum possible amplitude, depending on the frequency).

D. Participants

The experiment was performed with 28 participants. The participants ranged between 19 and 62 years, with a mean age of 28 years and a standard deviation of 11 years. Among the participants, 12 were women.

Participants were not compensated for their time. The study was approved by the research ethics committee of the University of Paris (IRB : 00012020-47).

III. RESULTS

A synthesis of the measured thresholds is shown in Figure 6 and Table II. The vibrotactile sensitivity seems to increase with frequency until it reaches a maximum at 160 Hz. Beyond this frequency, the sensitivity decreases constantly. Over these phases of decay and growth, the vibrotactile sensitivity seems to follow slopes of -6 and +9 dB per octave, respectively. The vibrotactile sensitivity curve is therefore U-shaped, centered at 160 Hz.

IV. DISCUSSION

It is a classical result of the literature that the vibrotactile sensitivity first increases in low frequencies, reaches a maximum, and then decreases in higher frequencies [START_REF] Verrillo | Investigation of some parameters of the cutaneous threshold for vibration[END_REF], [START_REF] Verrillo | Psychophysics of vibrotactile stimulation[END_REF], [START_REF] Verrillo | Sensation magnitude of vibrotactile stimuli[END_REF], [33], [START_REF] Brisben | Detection of vibration transmitted through an object grasped in the hand[END_REF]. The U-shape of this resulting curve (see Figure 6) is also classically observed [START_REF] Verrillo | Vibration sensation in humans[END_REF], [START_REF] Brisben | Detection of vibration transmitted through an object grasped in the hand[END_REF], [START_REF] Mills | Limitations of vibrotactile thresholds[END_REF], regardless of the stimulated body area. The maximum of sensitivity (minimum of the curve in Figure 6) is obtained here around 160 Hz, in agreement with thresholds obtained for vibrations parallel to the skin surface in this same body area [START_REF] Aevarsson | Vibrotactile threshold measurements at the wrist using parallel vibration actuators[END_REF]. It can therefore be assumed that the direction of application of the vibrations has little effect on the frequency of maximum sensitivity for this body area.

The inter-individual variability in the present experiment seems less important here than in other studies [START_REF] Aevarsson | Vibrotactile threshold measurements at the wrist using parallel vibration actuators[END_REF], [START_REF] Jansson | Vibrotactile thresholds on the mastoid and forehead position of deaf patients using radioear b71 and b81[END_REF], which may result from the difference in the direction of stimulation.

V. CONCLUSION

Perceptual thresholds in amplitude were estimated for vibrotactile signals applied normally to the internal side of the wrist. The trends observed on these perception thresholds (maximum sensitivity around 160 Hz, U-shaped curve) are in agreement with results already published, but measured with other vibrotactile devices and/or on a another body area.

The threshold search method used here is a classical method which allows a fine estimation of the threshold but leads to long experimental sessions. The estimation of a threshold for a frequency takes a few minutes. Some participants mentioned the monotony of the task. A fatigue phenomenon could have occurred and disturbed the results (the random order of the frequencies however probably limits its impact in the analysis). The use of a faster threshold search method (e.g. Audioscan [START_REF] Meyer-Bisch | Audioscan: A high-definition audiometry technique based on constant-level frequency sweeps -a new method with new hearing indicators[END_REF]) may limit this problem and allow thresholds to be obtained for a greater number of frequencies.

One part of the protocol could be improved to increase the reliability of the results by measuring the vibrator displacement for each participant. This measurement would allow to directly consider the actual displacement of the vibrator instead of estimating it from measurements made on a single person.

Now that the detection thresholds on the volar wrist have been estimated, a next step could be to estimate the isosensitivity curves. It was already shown on other body areas that they tend to remain parallel to the amplitude perception threshold curve [START_REF] Verrillo | Vibration sensation in humans[END_REF]. One can expect a similar trend for normal vibrations on the wrist, but this remains to be checked. Other psychophysical variables can also be studied in this context (JNDs in amplitude and frequency, for example). Variables that are fixed here (such as hand tightening and orientation, size of the surface in contact with the skin, type of transducers) may also be a focus of future work.

VI. APPENDIX : TIGHTENING FORCE ESTIMATION

Different methods of measuring the tightening force of the watchband were tested. The most classical one using a FSR sensor (Force Sensitive Resistor) was tested first. It is however impossible, in our case, to have a sufficient accuracy and fidelity of measurement. The forces involved in these tightening are too low to be measured repeatedly with a FSR sensor. Moreover, the introduction of this sensor between the vibrator and the skin modifies the coupling in an unknown way.

An original technique to indirectly estimate the force exerted by the tightening of the bracelet was designed. For this, the vibrator, without the watchband, is placed on the skin and different amounts of modeling clay are successively put on the vibrator to simulate different weights/tightening forces of the bracelet on the vibrator. The electrical impedance (Z = U/I where U is the input voltage and I the current) of the vibrator is measured for each added mass, the resonance frequency and the modulus of impedance at this frequency are extracted. The Thiele-Small model suggests that the addition of mass modifies these physical quantities [START_REF] Small | Passive-radiator loudspeaker systems part 1: Analysis[END_REF], [START_REF]Closed-box loudspeaker systems-part 2: Synthesis[END_REF].

A linear regression relating the value of the impedance modulus as a function of mass is determined (see dots and black line in Figure 7). From this regression, it is possible to estimate the force from the impedance modulus values. At the begining of the main experiment, the impedance modulus at resonance was measured, and the tightening estimated from linear regression (see gray squares in Figure 7). According to this method, the tightening of the watchband exerts a force on the vibrator of 1 N on average.

Fig. 1 .

 1 Fig. 1. Right-handed participant during the threshold measurement experiment. The vibrator is tightened by the wristband on the left volar wrist (corresponding to non preferred hand), the latter being placed on a foam limiting the transmission of vibrations. The participant uses their right hand to respond and identify the interval containing the stimulus.

Fig. 2 .

 2 Fig. 2. Vibrotactile wristband in place on the wrist of a participant and pointed by the laser vibrometer.

Fig. 3 .

 3 Fig. 3. Vibrator displacement as a function of frequency (different colors and markers) and input voltage (x-axis). The points represent averages over six measurements.

Fig. 4 .

 4 Fig. 4. Logigram of the threshold search method.

Fig. 5 .

 5 Fig. 5. Example of the threshold search method, for a given participant and frequency. Circles (resp. crosses) indicate correct (resp. wrong) answers. Black color indicates points separated by 3-dB steps, and gray color points separated by 1 dB steps.

Fig. 6 .

 6 Fig. 6. Threshold of perception for the 7 frequencies. The average of the thresholds is represented by a dot, the vertical lines represent an interval of 2 standard deviations around this average.

Frequency

TABLE I MAXIMUM

 I MEASURED THD, AS A FUNCTION OF THE SIGNAL'S FUNDAMENTAL FREQUENCY.

	Frequency (in Hz)	25	40	80	160	250	320	640
	Highest THD (in %)	14.8	10.2	8.3	5.8	2.5	1.8	1

TABLE II AMPLITUDE

 II DETECTION THRESHOLD AS A FUNCTION OF FREQUENCY.

DAEX19CT made by Dayton Audio (4Ω and 5W).

Internal face of the wrist.

Silicone watchband with magnetic closure by Tasikar.

https://uk.rs-online.com/web/p/ear-plugs/1905877

https://www.3m.co.uk/3M/en_GB/p/d/b00037368/

Scarlett 2i2 by Focusrite, 24 bits, 44.1 kHz.

FX1002A made by FX-AUDIO.

OFV-353 made by Polytec.

OFV-3001 made by Polytec.


Roberto Sabatini
email: sabatinr@erau.edu

Alessandro Monti

Yan Pailhas
email: yan.pailhas@cmre.nato.int

Angeliki Xenaki

Paul Cristini

An arbitrary-order immersed interface

This paper proposes an arbitrary-order immersed interface method for simulating the two-dimensional propagation of acoustic and elastic waves through fluid/solid interfaces. The present technique involves two main ingredients: 1) the linearized equations of continuum mechanics are simulated through an ADER (Arbitrary highorder schemes using DERivatives) scheme of arbitrary-order in both space and time (J. Comput. Phys 197(2): 532-539, 2004); 2) the jump conditions along the material interfaces are taken into account through the "explicit simplified interface method" (ESIM) derived by Lombard and Piraux (J. Comput. Phys 195(1): 90-116, 2004). In order to implement the ESIM, the computation of arbitrary-order spatial derivatives of the interface conditions is required. To this end, an algorithm for their numerical calculation, which does not require their explicit analytical expressions, is developed. Two numerical experiments involving flat and curved interfaces are finally discussed. When increasing the order of both the ADER scheme and of the interface treatment, the improvement of the convergence and of the accuracy of the numerical method is more specifically demonstrated by comparing the numerical results with analytical solutions.

I. INTRODUCTION

Simulating accurately the propagation of mechanical waves through heterogeneous media is essential in numerous scientific and engineering problems [START_REF] Finn | Computational Ocean Acoustics[END_REF] . In many physical applications, the propagation of acoustic and elastic waves is fairly well described by the linearized equations of continuum mechanics (LECM), i.e. the linearized Euler equations in fluids [START_REF] Bruneau | Equations of Motion in Dissipative Fluid[END_REF][START_REF] Pierce | The Wave Theory of Sound[END_REF] and the linear elastodynamic system in solids [START_REF] Ewing | Fundamental Equations and Solutions[END_REF] . Finite-difference (FD) methods are widely employed for solving such partial differential equations (PDEs) in single material domains. They are indeed straightforward to implement, inherently well-suited for parallelization techniques [START_REF] Sabatini | Three-dimensional direct numerical simulation of infrasound propagation in the Earth's atmosphere[END_REF] , and can reach very high accuracy [START_REF] Berland | High-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems[END_REF][START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF] . However, unlike finite-element methods [START_REF] Komatitsch | Wave propagation near a fluid-solid interface: A spectral-element approach[END_REF][START_REF] Zampolli | A computationally efficient finite element model with perfectly matched layers applied to scattering from axially symmetric objects[END_REF][START_REF] Zampolli | Low-to mid-frequency scattering from elastic objects on a sand sea floor: Simulation of frequency and aspect dependent structural echoes[END_REF][START_REF] Cristini | Some illustrative examples of the use of a spectralelement method in ocean acoustics[END_REF][START_REF]Influence of the sediment characteristics and of the level of burial on the acoustic response of a hollow cylinder in shallow water[END_REF][13] , FD schemes perform poorly in problems involving interfaces [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF] , such as the seafloor in underwater acoustics or the surface of an object immersed in a fluid. Firstly, FD formulas are valid only for sufficiently smooth functions, while, along a boundary between media, the mechanical fields are generally discontinuous or even defined differently in different adjacent domains. Secondly, the wave variables must verify not only the LECM but also appropriate constraints 2 that are not incorporated in FD algorithms. Finally, FD methods are typically implemented on regular Cartesian grids. Consequently, an interface is implicitly represented as a series of stairs and steps, which generates spurious reflected and transmitted waves and degrades the quality of the numerical solution unless sufficiently small grid spacings are employed [START_REF] Collino | Fictitious domain method for unsteady problems[END_REF] . To sum up, without an adequate treatment of the interfaces, the reflection and the transmission of acoustic and elastic waves are generally not correctly described by FD algorithms.

To cope with these issues, one of the strategies proposed in the literature employs the immersed interface methods (IIMs) [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF][START_REF] Zhang | The immersed interface method for acoustic wave equations with discontinuous coefficients[END_REF][START_REF] Farzi | High order immersed interface method for acoustic wave equation with discontinuous coeffcients[END_REF][START_REF] Piraux | A new interface method for hyperbolic problems with discontinuous coefficients: One-dimensional acoustic example[END_REF][START_REF] Lombard | Numerical modeling of elastic waves across imperfect contacts[END_REF][START_REF] Lombard | Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves[END_REF][START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF][START_REF] Chiavassa | Time domain numerical modeling of wave propagation in 2-D heterogeneous porous media[END_REF][START_REF] Abraham | A correction function method for the wave equation with interface jump conditions[END_REF][START_REF] Sabatini | An immersed interface method for the solution of the standard parabolic equation in range-dependent ocean environments[END_REF][START_REF] Singhal | A novel higher order compact-immersed interface approach for elliptic problems[END_REF][START_REF] Singhal | An efficient explicit jump high-order compact immersed interface approach for transient incompressible viscous flows[END_REF] . Among the IIMs, two particular algorithms for solving the linearized equations of continuum mechanics deserve consideration: the IIM developed by Zhang and LeVeque [START_REF] Zhang | The immersed interface method for acoustic wave equations with discontinuous coefficients[END_REF] and the method proposed by Lombard and Piraux [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF] . In the former, the coefficients of the FD formula are modified near the interfaces to account for the aforesaid constraints, thereby maintaining an acceptable level of accuracy. In the latter, a unique FD scheme is instead employed in the whole computational domain, and ghost or modified variables are used where the FD stencil crosses a boundary. Both techniques allow for the geometrical features of the interfaces. The latter method is however more attractive for two main reasons: firstly, it exhibits higher numerical stability, which is indispensable in configurations involving materials with high contrasts of the physical parameters; secondly, unlike the Zhang's IIM, it can possibly reach arbitrary orders of accuracy. Nevertheless, despite its potential, current implementations of the latter IIM are only second-or fourth-order accurate. Unfortunately, high-quality numerical simulations often require a high resolution in terms of points-per-wavelength, especially when the number of grid points in the computational domain is limited, for instance by the GPU memory.

The core of the algorithm proposed by Lombard and Piraux 14 lies in the computation of the ghost variables in the neighborhood of a boundary between media. The ghost variables are deduced from smooth extensions of the mechanical fields into the adjacent domains. This calculation is performed through the "explicit simplified interface method" (ESIM), which requires high-order spatial derivatives of the interface constraints. More specifically, Lombard and co-authors argued, on the basis of the work carried out by Gustafsson 27 , that conditions for the (k -1) th -order derivative of the unknown field are needed for a FD scheme to retain k th -order global accuracy [START_REF] Lombard | Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves[END_REF][START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF] . Unfortunately, deriving explicit expressions for such derivatives is an extremely tedious task, even for small values of k. For this reason, Lombard and Piraux 14 suggested the use of symbolic calculus software.

This manuscript aims at elaborating a formally arbitrary-order generalization of the twodimensional IIM proposed by Lombard and Piraux [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF] , where the main novelty is the automatic computation of arbitrary-order spatial derivatives of the interface conditions. Specifically, the present method consists in the coupling of two main ingredients:

1. an ADER (Arbitrary high-order schemes using DERivatives) scheme of order s in both space and time [START_REF] Toro | Towards Very High Order Godunov Schemes[END_REF][START_REF] Schwartzkopff | Fast high order ADER schemes for linear hyperbolic equations[END_REF][START_REF] Krause | Arbitrary high order finite volume schemes for linear wave propagation[END_REF] to compute the numerical solution of the linearized equations of continuum mechanics;

2. a k th -order ESIM for the treatment of the interfaces and the calculation of the ghost values in the neighborhood of material interfaces.

In order to simplify the forthcoming discussion, the proposed method is referred to as ADESIM sk. The main contribution of this work is the development of a novel algorithm that calculates the arbitrary-order spatial derivatives of the interface conditions. Besides, although the ADER schemes were originally derived in the framework of the finite-volume method [START_REF] Toro | Towards Very High Order Godunov Schemes[END_REF][START_REF] Schwartzkopff | Fast high order ADER schemes for linear hyperbolic equations[END_REF][START_REF] Krause | Arbitrary high order finite volume schemes for linear wave propagation[END_REF] , a derivation based on finite differences is additionally carried out [START_REF] Lombard | Modélisation numérique de la propagation et de la diffraction d'ondes mécaniques[END_REF] . Finally, the proposed algorithm would be beneficial for the extension of the method to three-dimensional configurations as it allows the explicit calculation by hand of a high number of terms for the spatial derivatives of the interface conditions to be avoided (for instance, see the appendix in Lombard [START_REF] Lombard | Méthodes numériques pour la propagation des ondes mécaniques et acoustiques en présence d'interfaces[END_REF] ).

The paper is organized as follows. The physical problem, the governing equations, and the interface conditions are introduced in Section II. Section III is devoted to the development of the novel algorithm for calculating arbitrary-order spatial derivatives of the interface conditions. The ADESIM sk method is then described in Section IV. In Section V, two numerical experiments are discussed to assess the validity of the method. The improvement of the numerical accuracy when increasing the order of both the numerical scheme and of the interface treatment is more particularly examined by comparing the numerical results with analytical solutions. Specifically, the numerical tests will show a significant improvement of the accuracy for high values of sk, e.g. s = 6 and k = 5, compared to the standard method with s = 2 and k = 2 [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF] . Concluding remarks are finally drawn in Section VI.

II. DESCRIPTION OF THE PHYSICAL PROBLEM

Introduce a Cartesian coordinate system Oxy and consider a generic interface Γ separating a fluid domain Ω 1 from a solid medium Ω 2 (cf. Figure 1). The curve Γ is described by a parameter τ as

Γ = {(x(τ ), y(τ )) | τ ∈ I ⊆ R} , (1)
and its tangential and normal vectors, labeled as t and n, are defined by the expressions

t(τ ) = dx dτ dy dτ T n(τ ) = dy dτ - dx dτ T ( 2
)
where the symbol T indicates the matrix transposition. The interface Γ is of class C k , i.e the functions d k x/dτ k and d k y/dτ k are everywhere continuous up to a order k.

In the study reported in this paper, the mechanical waves propagating through the media Ω 1 and Ω 2 are considered as small-amplitude perturbations of initial undisturbed ambient states (linearity assumption). The fluid ambient state is characterized by the density ρf and by the speed of sound cf . The acoustic propagation is described by the two components of the velocity perturbation u, u f

x and u f y , and by the pressure perturbation p, collected in the vector

U = u f x u f y p T .
The physical parameters for the solid ambient state are the density ρs , the speed of the P -waves cp and the speed of the S-waves cs , both expressed in terms of the Lamé coefficients, λ and μ33 , i.e. cp = ( λ + 2μ)/ρ s and cs = μ/ρ s . The elastic propagation is described by the two components of the elastic velocity u, u s x and u s y , and by the three components of the elastic stress tensor Σ, σ xx , σ xy and σ yy , collected in the vector U = u s x u s y σ xx σ xy σ yy T .

Along with the small perturbation hypothesis, the following premises are made:

• within each material domain, all ambient variables are independent of space and time;

• the fluid is at rest;

• the solid is isotropic and elastic;

• attenuation phenomena and body-force effects (such as the effect of gravity) are negligible.

Under these assumptions, the linearization of the equations of continuum mechanics leads to a first-order linear hyperbolic system [START_REF] Pierce | The Wave Theory of Sound[END_REF][START_REF] Ewing | Fundamental Equations and Solutions[END_REF][START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF] :

∂U ∂t = -A ∂U ∂x -B ∂U ∂y . (3)
In Eq. ( 3), the matrices A and B are constant and depend on the physical parameters. The matrices are defined in a space A, B ∈ R 3×3 in the fluid, and A, B ∈ R 5×5 in the solid. Finally, for the mathematical problem to be well-posed, the following conditions must be additionally verified at a generic point P on the interface Γ and at each instant of time • the normal velocity is continuous,

u T n = 0, ∀P ∈ Γ, ∀t; (4a)
• the normal stress is continuous,

n T Σn = -p, ∀P ∈ Γ, ∀t; (4b)
• the tangential stress vanishes,

t T Σn = 0, ∀P ∈ Γ, ∀t. (4c)

III. ARBITRARY HIGH-ORDER INTERFACE CONSTRAINTS

In this section, the algorithm developed to compute arbitrary-order spatial derivatives of the interface conditions (4) is described. Such derivatives are employed in the ADESIM sk method presented in Sec. IV.

A. Notations

Consider a point P on Γ (cf. Figure 1). The limits of the solution U (x, y, t) and of its spatial derivatives up to the order k are collected in the vector

U k i = lim M →P,M ∈Ω i U , ∂U ∂x , ∂U ∂y , ..., ∂ α U ∂x α-β ∂y β , ..., ∂ k U ∂y k T , i = 1, 2, α, β = 0, ..., k, (5)
where the subscript i indicates the medium. By definition, the n th term in U k i is

U k i [n] = lim M →P,M ∈Ω i ∂ α U [l] ∂x α-β ∂y β , (6)
where U [l] is the l th component of U . The index n is identified by the function

n = idx(α, β, l, d) = dα(α + 1) 2 + dβ + l, l ∈ 1, . . . , d, (7)
where d is the length of U . The total number of terms in

U k i is simply idx(k, k, d, d) = dk(k + 1) 2 + dk + d. (8)
As a result, the vectors U k 1 and U k 2 contain respectively 3(k 2 +3k +2)/2 and 5(k 2 +5k +2)/2 components.

For the sake of consistency, the following notations are employed for the matrices A and B:

A i = lim M →P,M ∈Ω i A, B i = lim M →P,M ∈Ω i B. (9)

B. Arbitrary-order derivatives of the interface constraints

The interface constraints (4) are rewritten as

   C 0 1 U 0 1 = C 0 2 U 0 2 L 0 2 U 0 2 = 0 , (10)
where the matrices C 0 1 , C 0 2 and L 0 2 are defined by the following expressions:

C 0 1 =      - dy dτ dx dτ 0 0 0 - dx dτ 2 - dy dτ 2      C 0 2 =      - dy dτ dx dτ 0 0 0 0 0 dy dτ 2 -2 dx dτ dy dτ dy dτ 2      L 0 2 = 0 0 dx dτ dy dτ dy dτ 2 - dx dτ 2 - dx dτ dy dτ . ( 11
)
Arbitrary-order interface conditions are developed by deriving Eqs. (10) with respect to the variables τ and t. The term

∂ k T 0 i U 0 i /∂τ k-q ∂t q , with T 0 i equal to either C 0 1 , C 0 2 , or L 0 2
, is given by the formula (for instance, see Hardy [START_REF] Hardy | Combinatorics of partial derivatives[END_REF] and references therein):

∂ k T 0 i U 0 i ∂τ k-q ∂t q = k-q q 1 =0 q q 2 =0 k -q q 1 q q 2 ∂ q 1 +q 2 T 0 i ∂τ q 1 ∂t q 2 ∂ k-q 1 -q 2 U 0 i ∂τ k-q-q 1 ∂t q-q 2 . ( 12
)
Since the matrix T 0 i only depends on the parameter τ , the double sum (12) does not vanish only if q 2 = 0. Therefore, Eq. ( 12) is simplified as

∂ k T 0 i U 0 i ∂τ k-q ∂t q = k-q l=0 k -q l d l T 0 i dτ l ∂ k-l U 0 i ∂τ k-q-l ∂t q . ( 13
)
The derivatives d l T 0 i /dτ l are computed through the formulas for the derivatives of products provided by Hardy [START_REF] Hardy | Combinatorics of partial derivatives[END_REF] . The function ∂ q U 0 i /∂t q in (13) is calculated through the Lax-Wendroff (or Cauchy-Kovalevskaya) procedure [START_REF] Schwartzkopff | Fast high order ADER schemes for linear hyperbolic equations[END_REF] , which consists in replacing, in Eqs. (3), the q th -order temporal derivative with spatial derivatives:

∂ q U 0 i ∂t q = (-1) q A i ∂ ∂x + B i ∂ ∂y q U 0 i . (14)
By employing combinatorial rules for the binomial term, the above system is further developed as

∂ q U 0 i ∂t q = q m=0 (-1) q M i,m,q ∂ q U 0 i ∂x q-m ∂y m , (15)
where M i,m,q is the sum of the n up unique permutations of (qm) matrices A i and m matrices B i . As an example, M i,1,3 is equal to the sum of the unique permutations of 2 matrices A i and

1 matrix B i , i.e. A i A i B i , A i B i A i , and B i A i A i , i.e. M i,1,3 = A i A i B i + A i B i A i + B i A i A i .
Inserting Eq. ( 15) in Eq. ( 13) yields the following expression:

∂ k T 0 i U 0 i ∂τ k-q ∂t q = k-q l=0 q m=0 (-1) q k -q l d l T 0 i dτ l M i,m,q ∂ q ∂x q-m ∂y m ∂ k-q-l U 0 i ∂τ k-q-l . ( 16
)
To calculate the term ∂ k-q-l U /∂τ k-q-l , a short digression is needed. Suppose that f, g : R → R are functions admitting n derivatives. According to Mennucci [START_REF] Mennucci | An intuitive presentation of faà di bruno's formula[END_REF] , the N th -derivative of the composition f (g(τ )) is given by the Faà di Bruno's formula 36 :

d N f (g(τ )) dτ N = N n=1 1 n! d n f (g) dg n (s 1 ,...,sn)∈C N N s 1 , s 2 , . . . , s n n r=1 d sr g(τ ) dτ sr . (17)
In Eq. ( 17), the internal sum is over all the compositions C N of the integer N into n parts, i.e. over all the n-tuples of positive indices s r , r = 1, . . . , n, that satisfy the constraint

s 1 + • • • + s n = N .
For given n and N , the total number of possible compositions is

N n=1 N -1 n -1 = 2 N -1 . (18)
By moving the derivative d n f (g)/dg n inside the product sign, Expression ( 17) is rewritten as

d N f (g(τ )) dτ N = N n=1 1 n! (s 1 ,...,sn)∈C N N s 1 , s 2 , . . . , s n n r=1 d sr g(τ ) dτ sr d dg f. (19)
According to Mishkov 37 , when a composite function f (g(τ )) depends on a vector variable g, g(τ ) = [g 1 (τ ), g 2 (τ ), . . . ], the derivative d sr g(τ )/dτ sr d/dg must be replaced by the operator (d sr g/dτ sr • ∇). As a consequence, the function ∂ k-q-l U i /∂τ k-q-l is equal to

∂ k-q-l U 0 i ∂τ k-q-l = k-q-l n=1 1 n! (s 1 ,...,sn)∈C k-q-l k -q -l s 1 , s 2 , . . . , s n n r=1 d sr x dτ sr ∂ ∂x + d sr y dτ sr ∂ ∂y U 0 i . ( 20
)
The derivative ∂ k T 0 i U 0 i /∂τ k-q ∂t q is then given by the formula

∂ k T 0 i U 0 i ∂τ k-q ∂t q = k-q l=0 q m=0 k-q-l n=1 (s 1 ,...,sn)∈C k-q-l (-1) q n! k -q l k -q -l s 1 , . . . , s n × × d l T 0 i dτ l M i,m,q ∂ q ∂x q-m ∂y m n r=1 d sr x dτ sr ∂ ∂x + d sr y dτ sr ∂ ∂y U 0 i . (21)
In Eq. ( 21), the factor between brackets is the product of n binomials of the form (a r + b r ), with r = 1, . . . , n. By virtue of the rule of product, this factor is rewritten as the sum of 2 n terms, each of which is the product of n monomials chosen in the n distinct aforesaid binomials.

Despite its apparent complexity, Eq. ( 21) is recast in the form

T q,k i U k i = ∂ k T 0 i U 0 i ∂τ k-q ∂t q . ( 22
)
Combining all the derivatives up to the order k yields the expression:

T 0 i U 0 i ∂T 0 i U 0 i ∂τ ∂T 0 i U 0 i ∂t . . . ∂ k T 0 i U 0 i ∂t k T = T 0,0 i T 0,1 i T 1,1 i . . . T k,k i T U k i = T k i U k i . ( 23
)
Enforcing the continuity of the derivatives of C 0 i U 0 i and requiring that the derivatives of L 0 2 U 0 2 vanish on Γ finally leads to the following system:

C k 1 U k 1 0 = 0C k 2 U k 2 0; (24a) L k 2 U k 2 0 = 00C k 2 U k 2 . ( 24b
)
For a given value of k, there are (k + 1) possible mixed derivatives. Consequently, the matrices C k i and L k i have respectively (k 2 + 3k + 2)/2 × 2 and (k 2 + 3k + 2)/2 × 1 rows, and

C k i ∈ R (k 2 +3k+2) × idx(k,k,d,d) , L k i ∈ R (k 2 +3k+2)/2 × idx(k,k,d,d) .

C. Irrotationality of the fluid flow and compatibility conditions for the elastic stresses

As highlighted by Lombard and Piraux 14 , the components of U are not independent, but linked by compatibility relationships. Such constraints, along with their higher-order derivatives, are employed as supplementary constraints for the treatment of the interfaces. For the sake of completeness, they are briefly recalled in the forthcoming lines.

Under the assumptions made in Section II, the fluid flow is irrotational. Therefore, the curl of the acoustic velocity vector must be identically nil:

∂u f x ∂y - ∂u f y ∂x = 0. ( 25
)
Deriving Expression (25) (kq -1) times with respect to x and q times with respect to y yields the following equation:

∂ k u f x ∂x k-q-1 ∂y q+1 - ∂ k u f y ∂x k-q ∂y q = 0, k ≥ 1, 0 ≤ q ≤ (k -1). (26)
For a given order k, (k 2 + k)/2 constraints are derived. Moreover, they must be verified at each location (x, y) in the fluid and in particular at the interface point P . Consequently, the vector U k 1 must satisfy a linear system of the form

F k 1 U k 1 = 0. (27)
The matrix

F k 1 ∈ R (k 2 +k)/2×idx(k,k,3,3
) is computed through the algorithm provided in Algorithm 1.

In a two-dimensional isotropic elastic solid, the three components of the stress tensor are completely determined by the two displacements along the horizontal and vertical axes. Hence, the functions σ xx , σ xy and σ yy are not independent and must verify a compatibility condition. More specifically,

α 2 ∂ 2 σ xx ∂x 2 + α 1 ∂ 2 σ yy ∂x 2 - ∂ 2 σ xy ∂x∂y + α 1 ∂ 2 σ xx ∂y 2 + α 2 ∂ 2 σ yy ∂y 2 = 0. (28)
The coefficients α 1 and α 2 are given in terms of the Lamé parameters, i.e. α 1 = ( λ + 2μ)/[4( λ + μ)] and α 2 = -λ/[4( λ + μ)]. Deriving Expression (28) (kq -2) times with respect to x and q times with respect to y yields the equation

α 2 ∂ k σ xx ∂x k-q ∂y q + α 1 ∂ k σ yy ∂x k-q ∂y q - ∂ k σ xy ∂x k-q-1 ∂y q+1 + α 1 ∂ k σ xx ∂x k-q-2 ∂y q+2 + α 2 ∂ k σ yy ∂x k-q-2 ∂y q+1 = 0, k ≥ 2, 0 ≤ q ≤ (k -2). (29)
For a given order k, (k 2k)/2 constraints can then be derived. As in the case of the fluid domain, they must be verified at each location (x, y) in the solid and in particular at the interface point P . As a consequence, the vector U k 2 must satisfy a linear system of the form

F k 2 U k 2 = 0. (30)
The matrix

F k 2 ∈ R (k 2 -k)/2×idx(k,k,5,5
) is calculated through the algorithms provided in Algorithm 2.

1: q = 1; 2: for n = 1 : k do 3:

for i = 0 : (n -1) do 4: F k 1 (q, idx(n, i + 1, 1, 3)) = +1; 5:
F k 1 (q, idx(n, i + 0, 2, 3)) = -1;

6:

q = q + 1;

7:
end for 8: end for Algorithm 1: Algorithm for the computation of the matrix F k 1 .

1: q = 1; 2: for n = 2 : k do 3:

for i = 0 : (n -2) do 4:

F k 2 (q, idx(n, i + 0, 3, 5)) = +α 2 ;

5:

F k 2 (q, idx(n, i + 0, 5, 5)) = +α 1 ;

6:

F k 2 (q, idx(n, i + 1, 4, 5)) = -1;

7:

F k 2 (q, idx(n, i + 2, 3, 5)) = +α 1 ;

8:

F k 2 (q, idx(n, i + 2, 5, 5)) = +α 2 ;

9:

q = q + 1;

10:

end for 11: end for Algorithm 2: Algorithm for the computation of the matrix F k 2 .

The constraints ( 27) and (30) are then used to reduce the number of unknowns in U k 1 and U k 2 , respectively. Let U k i,K and U k i,R be the vectors containing the degrees of freedom to Keep and Remove from U k i , respectively. In the fluid vector, the degrees of freedom ∂ n u y /∂x n-q ∂y q , for n = 1, ..., k and q = 0, ..., (n -1), are removed. For the solution in the solid medium, the unknowns ∂ n σ yy /∂x n-q ∂y q , for n = 2, ..., k and q = 0, ..., (n -2), are considered dependent variables. As a consequence, the vectors U k 1,K and U k 2,K have (k 2 + 4k + 3) and (2k 2 + 8k + 5) components, respectively. The system F k i U k i = 0 is conveniently rearranged as

F k i,K F k i,R   U k i,K U k i,R   = 0. (31)
Therefore, the relation

U k i,R = -F k i,R -1 F k i,K U k i,K holds and the vector U k i is conveniently rewritten as U k i = G k i U k i,K , (32)
where the matrices

G k 1 ∈ R idx(k,k,3,3) × (k 2 +4k+3) and G k 2 ∈ R idx(k,k,5,5) × (2k 2 +8k+5
) are defined by the formulas

G k 1 =   I k 2 +4k+3 -F k 1,R -1 F k 1,K   , G k 2 =   I 2k 2 +8k+5 -F k 2,R -1 F k 2,K   , (33)
where I q ∈ R q × q is the identity matrix of size q × q. Inserting Expression (32) in Eq. (24b) yields

L k i G k i U k i,K = 0, (34)
where the relation (24b) for i = 1 is verified by default being a fluid at rest. Therefore, the vector

U k i,K belongs to the kernel Z L k i G k i of L k i G k i . Let Z L k i G k
i be the matrices whose columns are the basis vectors of

Z L k i G k i . Then, any U k i,K ∈ Z L k i G k i can be expressed as U k i,K = Z L k i G k i V k i,K . (35)
The vectors V k 1 and V k 2 have respectively (k 2 + 4k + 3) and (3k 2 + 13k + 8)/2 components, and, consequently,

Z L k 1 G k 1 ∈ R (k 2 +4k+3) × (k 2 +4k+3) and Z L k 2 G k 2 ∈ R (2k 2 +8k+5) × (3k 2 +13k+8)/2 . The matrix Z L k i G k i is deduced from a Singular Value Decomposition (SVD) of L k i G k i . The interface conditions C k 1 U k 1 = C k 2 U k 2 are finally rewritten as S k 1 V k 1 = S k 2 V k 2 , (36)
where the matrices S k 1 and S k 2 are defined by the following formulas:

S k 1 = C k 1 G k 1 Z L k 1 G k 1 ∈ R (k 2 +3k+2) × (k 2 +4k+3) ; S k 2 = C k 2 G k 2 Z L k 2 G k 2 ∈ R (k 2 +3k+2) × (3k 2 +13k+8)/2 . (37)
To employ the ESIM, the solution V k 2 must be expressed as a function of V k 1 , and viceversa. System (36) is underdetermined and is solved through the SVD method. Its general solution is

V k 2 = S k 2 + S k 1 V k 1 + I (3k 2 +13k+8)/2 -S k 2 + S k 2 Λ k 1 , (38)
where S k 2 + is the pseudoinverse of S k 2 and Λ k 1 ∈ R (3k 2 +13k+8)/2×1 is any real column vector. An analogous result holds for V k 1 ,

V k 1 = S k 1 + S k 2 V k 2 + I (k 2 +4k+3) -S k 1 + S k 1 Λ k 2 , (39)
where

S k 1 + ∈ R (k 2 +4k+3) × (k 2 +3k+2) is the pseudoinverse of S k 1 and Λ k 2 ∈ R (k 2 +4k+3
)×1 is any real column vector.

IV. ADESIM sk METHOD

In this section, the ADESIM sk method is presented. The ADER scheme employed for solving the governing system (3) is first outlined. The ESIM procedure, with the novel algorithm for computing arbitrary-order spatial derivatives of the interface conditions outlined in Section III C, is then recalled for the sake of completeness.

A. ADER schemes for the governing system (3) Define a time step ∆t and a uniform Cartesian grid, consisting of N x and N y nodes along the xand y-axes, respectively (cf. Figure 2). The spacings are labeled as ∆x and ∆y and, for the sake of simplicity, they are assumed equal throughout this paper, ∆x = ∆y. The interface Γ is immersed in the computational mesh. Consider a generic point M of coordinates x i = i∆x and y j = j∆y, with i = 1, . . . , N x and j = 1, . . . , N y , and suppose temporarily that M is sufficiently far from Γ (blue node in Figure 2). The function U (x, y, t) is therefore smooth in a neighborhood of M . The term U (x i , y j , t n+1 ), evaluated at the instant t n+1 = (n + 1)∆t, can be formally computed through a Taylor series around t n :

U (x i , y j , t n+1 ) = U (x i , y j , t n ) + ∞ q=1 ∆t q q! ∂ q U ∂t q x=x i ,y=y j ,t=tn . (40)
Using Eq. ( 3) to replace the temporal derivative by spatial derivatives and employing the combinatorial rules for the binomial term, as done in Eq. ( 15), yields the expression

U (x i , y j , t n+1 ) = U (x i , y j , t n ) + s q=1 q m=0
(-1) q ∆t q q! M m,q ∂ q U ∂x q-m ∂y m x=x i ,y=y j ,t=tn + O(∆t s+1 ).

(41) Note that the temporal Taylor series has been truncated at the order s, with only even s considered in this work. The matrix M m,q is defined as in Section III. A numerical scheme is then obtained by replacing the spatial derivatives with finite-difference formulas. To this end, a square stencil centered at the point (x i , y j ) and employing at most (s + 1) × (s + 1) nodes is used. The approximation U n+1 i,j of the solution U (x i , y j , t n+1 ) is written as

U n+1 i,j = U n i,j + s/2 α=-s/2 s/2 β=-s/2 D α,β U n i+α,j+β . (42)
The numerical scheme (42) has the structure of a single-step finite difference method, with a stencil of (s + 1) × (s + 1) nodes centered around the point M , and can be abstractly rewritten as

U n+1 i,j = H U n i-s/2,j-s/2 , . . . , U n i+s/2,j+s/2 , (43)
* O y x M I (x I , y J ) M Q P Q (x i , y j ) M (x i , y j ) ∆x ∆y

FIG. 2. (color online) Computational grid

where H is a linear discrete operator. The scheme H is of order s in space and in time, i.e. the local truncation error behaves as O(∆t s+1 , ∆x s+1 , ∆y s+1 ). Moreover, the grid spacing ∆x is related to the time step ∆t by the Courant-Friedrichs-Lewy number

CFL = max(c f , cp , cs ) ∆t ∆x , (44)
whose maximum value is 1/ √ 2 for s = 2 and 1 for s ≥ 4 [START_REF] Krause | Arbitrary high order finite volume schemes for linear wave propagation[END_REF] .

B. Explicit simplified interface method (ESIM)

Consider at present a node M I of coordinates (x I , y J ) and close to the interface Γ (cf. Figure 2). Assume that one of the stencil points Q of coordinates (x I+α , y J+β ), α , β ∈ [-s/2, s/2], α = 0, β = 0, lies in the adjacent domain (red asterisk in Figure 2). The node M I is therefore regarded as irregular and Formula (43) cannot be applied, without modifications, to calculate the approximation U n+1 I,J of U (x I , y J , t n+1 ). The ESIM method provides a strategy to overcome this issue. Its basic idea is to employ Expression (43) at the node M I as well, upon replacing the term U n I+α ,J+β with an appropriate modified or ghost value U n I+α ,J+β . As an illustration, the numerical solution U n+1 I,J at the red node displayed in Figure 2 would be computed as

U n+1 I,J = H (U n I-1,J-1 , U n I+0,J-1 , U n I+1,J-1 , H (U n I-1,J+0 , U n I+0,J+0 , U n I+1,J+0 ,
H (U n I-1,J+1 , U n I+0,J+1 , U n I+1,J+1 ) . ( 45
)
Let Ω i and Ω i be the domain of the irregular node M I and the opposite medium, respectively. Moreover, let P Q be the orthogonal projection of Q on the interface Γ, with coordinates (x P Q , y P Q ). The ghost value U n I+α ,J+β can be interpreted as the point value of a smooth extension U of the mechanical field U in Ω i into the domain Ω i (the fluid solution is smoothly extended in the solid side and viceversa). Such an extension is defined by a k th -order Taylor series around the point P Q seen from the domain Ω i :

U (x, y, t) = Π k i (x, y)U k i (x P Q , y P Q , t). ( 46
)
In Eq. ( 46), the components of the matrix d) are given by the expression

Π k i ∈ R d × idx(k,k,d,
Π k i [l, idx(v, w, l, d)](x, y) = (x -x P Q ) v-w (y -y P Q ) w (v -w)!w! , l = 1, . . . , d, r, s = 1, . . . , k. ( 47
)
The modified value U n I+α ,J+β is then calculated as

U n I+α ,J+β = Π k i (x I+α , y J+β )U k,n i,P Q , ( 48
)
where U k,n i,P Q is an approximation of the vector

U k i (x P Q , y P Q , t n ). Therefore, computing U n
I+α ,J+β is tantamount to estimate U k,n i,P Q . Assume at present that the point Q is located in the solid side, as in Figure 2. In order to evaluate U k,n i,P Q , a set B of grid nodes surrounding P Q is employed (cf. Figure 3). This

O y x Q P Q Γ B B 1 B 2 Ω 2 Ω 1 FIG. 3. ( color
online) Set of nodes B for the estimation of the vector U k,n i,P Q set is enclosed in a circle centered at the point P Q with a radius r B and is divided into two subsets B i , i = 1, 2, by the interface Γ. The solution at a node of indexes (i, j) in the subset B 1 is given, at the instant t n , by the expression

U (x i , y j , t n ) = Π k 1 (x i , y j )U k 1 (t n ) + O(∆x k+1 ) = Π k 1 (x i , y j )G k 1 Z L k 1 G k 1 V k 1 (t n ) + O(∆x k+1 ), (49)
which is rewritten as

U (x i , y j , t n ) = Π k 1 (x i , y j )G k 1 Z L k 1 G k 1 I (k 2 +4k+3) O Λ k 1 V k 1 (t n ) Λ k 1 (t n ) + O(∆x k+1 ), (50)
where O Λ k 1 is the null matrix of size (k 2 + 4k + 3) × (3k 2 + 13k + 8)/2. At a point of indexes (i, j) located in the subset B 2 , the solution at the instant t n is instead given by the following formula:

U (x i , y j , t n ) = Π k 2 (x i , y j )U k 2 (t n ) + O(∆x k+1 ) = Π k 2 (x i , y j )G k 2 Z L k 2 G k 2 V k 2 (t n ) + O(∆x k+1 ). (51)
Using Eq. ( 38) to recast the vector V k 2 as a function of the vector V k 1 yields

U (x i , y j , t n ) = Π k 2 (x i , y j )G k 2 Z L k 2 G k 2 S k 2 + S k 1 (I (3k 2 +13k+8)/2 -S k 2 + S k 2 ) V k 1 (t n ) Λ k 1 (t n ) +O(∆x k+1 ).
(52) Combining Equation (50) and Equation ( 52) for all the nodes in B leads to an expression of the following form:

(U n ) B ≡     U (x i 1 , y j 1 , t n ) U (x i 2 , y j 2 , t n ) . . .     = M B,1 (x i 1 , y j 1 , x i 2 , y j 2 , . . . )     V k 1 (t n ) Λ k 1 (t n )     +     O(∆x k+1 ) . . . O(∆x k+1 )     . ( 53
)
The radius r B of B is chosen in such a way that System (53) is overdetermined. A convenient restriction M + B,1 of the pseudoinverse M B,1 of the matrix M B,1 provides an approximation of the unknown vector V k 1 and

U k,n i,P Q = G k 1 Z L k 1 G k 1 M + B,1 (U n ) B . (54)
The modified value U n I+α ,J+β is finally calculated as

U n I+α ,J+β = Π k i (x I+α , y J+β )G k 1 Z L k 1 G k 1 M + B,1 (U n ) B . (55)
For a node Q in the fluid side, an analogous procedure leads to a formula of the form

U n I+α ,J+β = Π k i (x I+α , y J+β )G k 2 Z L k 2 G k 2 M + B,2 (U n ) B . (56)
The computations of U k,n i,P Q and U n I+α ,J+β involve the high-order constraints described in Section III. As a result, the modified or ghost value U n I+α ,J+β implicitly allows taking into account the interface conditions, which are necessary for a correct description of the reflection and transmission phenomena of waves at a boundary between different media. The above procedure must be repeated at each time step for all irregular nodes and grid points located in the opposite side of Γ, before using the ghost values in the scheme (43).

V. NUMERICAL EXPERIMENTS

Two numerical experiments are carried out. In the first test case, the propagation of a plane wave through a flat fluid/solid interface is simulated. This benchmark aims at demonstrating the improvement of the accuracy and of the convergence rate when increasing the orders s and k of the ADER scheme and the interface treatment. The second experiment concerns the acoustic scattering of a cylindrical pulse by a water-filled aluminum cylindrical shell immersed in water. Its goals are to show the ability of the ADESIM sk method to handle curved interfaces and to examine the resolutions in terms of points-per-wavelength for various values of s and k.

A. First test case

In the first test case, the material boundary is described by a parameter τ as

Γ = (x, y) ∈ R 2 | x(τ ) = 0, y(τ ) = τ, τ ∈ R , (57)
and the physical properties of the media are set to ρf = 1000 kg m -3 , cf = 1500 m s -1 , ρs = 2600 kg m -3 , cp = 4000 m s -1 , and cs = 2000 m s -1 . The initial field is a plane wave propagating in the fluid with an angle ϑ of 21 • with respect to the horizontal axis and is defined by the expression

U (x, y, t) = - cos(ϑ) cf sin(ϑ) cf ρf T F t - x c 1 , (58)
for t = 0. The function F is a C 6 spatially-bounded combination of truncated sinusoids and is given by the formula 21

F(ξ) =      4 m=1 a m sin(b m ω 0 ξ) ξ ∈ 0, 2π ω c 0 otherwise , (59)
where a 1 = 1, a 2 = -21/32, a 3 = 63/768, a 4 = -1/512, b m = 2 m-1 , and ω 0 = 2π × 8 kHz. The physical domain of interest spans the square [-2, 2] m × [-2, 2] m. At its boundaries, the analytical solution provided by Lombard 32 is employed to update the mechanical field U . Different simulations are performed with various values of ∆x, s, and k. The Courant-Friedrichs-Lewy number CFL is set equal to 0.95 for s ≥ 4 and to 0.95/1.5 for s = 2. The numerical results are compared with the analytical solution.

The stresses σ yy obtained at three different instants of time with ∆x = 1/256 m, s = 6 and k = 5 are reported in Figure 4 (in the fluid σ yy = -p). In the fluid, the mechanical field consists of an incident and of a reflected acoustic waves. In the solid, transmitted P-and S-waves are observed.

The errors in norm L ∞ calculated at the instant t = 7.4 × 10 -4 s, on the line y = -0.1 m, are displayed in Figure 5 with the convergence rates. In agreement with the theoretical results of Gustafsson 27 , a (s -1) th -order interface treatment coupled with a s th -order ADER scheme is sufficient to maintain s th -order global accuracy. The schemes ADESIM x (m)

-2 -1 +0 +1 +2 y (m) Ω 1 Ω 2 Γ (a) t = 0.37 ms -2 -1 +0 +1 +2
x (m) x (m) improves the numerical accuracy.

B. Second test case

In the second test case, the internal and external fluid/solid interfaces are described by the parametric equations

Γ internal = {(x, y) ∈ R 2 | x(τ ) = r i cos(τ ) + 0.5, y(τ ) = r i sin(τ ), τ ∈ [0, 2π[} , Γ external = {(x, y) ∈ R 2 | x(τ ) = r e cos(τ ) + 0.5, y(τ ) = r e sin(τ ), τ ∈ [0, 2π[} , (60)
where r i = 0.4 m and r e = 0.5 m. The density and the speed of sound of the internal and external fluids are ρf = 1000 kg m -3 and cf = 1500 m s -1 . The properties of the solid are ρs = 2700 kg m -3 , cp = 6420 m s -1 , and cs = 3040 m s -1 . The initial mechanical field is nil and an incident cylindrical pulse (Ricker wavelet) is excited in the external fluid by the pressure source

S(x, y, t) = - 1 πα 2 1 - ω 2 c 2 (t -t c ) 2 e -ω 2 c (t-tc) 2 /4 e -(x-xc) 2 /α 2 -(y-yc) 2 /α 2 [0 0 1] T , (61)
where ω c = 2π × 20 kHz, t c = 0.1 ms, x c = -0.5 m, y c = 0 m, and α = 0.005 m. The term (61) is integrated in the numerical solution at the end of each time step through the Euler method. The computational frame spans the square [-1, 2] m × [-1.5, 1.5] m. At its boundaries, the radiation conditions proposed by Tam and Webb 38 are employed to minimize the reflections of outgoing waves. Simulations are performed with ∆x = 3.75 × 10 -3 m by employing the following schemes: ADESIM 2 -2, ADESIM 4 -3, ADESIM 6 -3, and ADESIM 6 -5. The Courant-Friedrichs-Lewy number CFL is set equal to 0.95 for s ≥ 4 and to 0.95/1.5 for s = 2. The numerical results are compared with the analytical solution. The latter is computed via standard techniques for acoustic scattering by cylindrical shells [START_REF] Alexander | Scattering of low-frequency sound by compact objects in underwater waveguides[END_REF] , with the incident acoustic field calculated as suggested in Sabatini et al. [START_REF] Sabatini | Far-field analytical solutions of the non-homogeneous helmholtz and wave equations for spatially non-localized sources[END_REF] .

The stress fields σ xx obtained at four different instants of time through the scheme ADESIM 6 -5 are displayed in Figure 7. At t = 0.33 ms (cf. Figure 7(a)), a cylindrical impulsive wave has been excited in the ambient fluid and is travelling toward the shell. At t = 0.66 ms (cf. Figure 7(b)), the pulse has been reflected by the external interface and is propagating through the object. At t = 1.32 ms (cf. Figure 7(c)) and t = 2.74 ms (cf. Figure 7(d)), the scattered field has filled the whole physical domain of interest.

The signals of pressure calculated at the point of coordinates (-0.5, 0) m are shown in Figures 8(a) and 8(b) along with the analytical solution. The corresponding one-sided energy spectral densities are displayed in Figure 8(c). While numerical dispersion is observed on the results computed with the second-and the fourth-order ADESIM schemes, the curves calculated through the sixth-order methods are in excellent agreement with the exact solution. This trend is also visible in the ESD diagram, which clearly shows the increase of the numerical accuracy with the orders s. As an illustration, for s = 6, appreciable differences between the numerical and the exact ESDs only appear for frequencies higher than about 52 kHz. In order to quantify the accuracies of the computations, a relative error E FT{-p} for the modulus of the Fourier transform of a numerical time series is introduced:

E FT{-p} (f ) = |FT{-p numerical }(f ) -FT{-p analytical }(f )| |FT{-p analytical }(f )| . ( 62
)
For each numerical signal, the function E FT{-p} (f ) is smoothed through a moving average with a window length of 2 kHz. The resulting curves are plotted in Figure 8(d). For each scheme, the smallest resolved wavelength λ min is then defined as the minimum value of λ = cf /f such that E FT{-p} (f ) < 10 -1 . The results are reported in Table I. Note that with ADESIM 6 -5, the largest resolved frequency f max is 342% higher than the one resolved with standard method ADESIM 2 -2 14 . Consequently, the number of points required for resolving the smallest wavelength λ min drastically reduces by the same percentage amount, bringing to a significant gain in computational memory requirements.

Increasing the orders of both the ADER scheme and the interface treatment, while keeping the grid spacing ∆x constant, clearly enhances the numerical resolution. As an illustration, a sixth-order ADER method coupled with a fifth-order ESIM achieves a resolution of about λ min /∆x = 7.6 points-per-wavelength, which is comparable to the accuracy of very-highquality solvers, such as the well-known SPECFEM code based on the spectral element method [START_REF] Komatitsch | Wave propagation near a fluid-solid interface: A spectral-element approach[END_REF][START_REF] Cristini | Some illustrative examples of the use of a spectralelement method in ocean acoustics[END_REF][START_REF]Influence of the sediment characteristics and of the level of burial on the acoustic response of a hollow cylinder in shallow water[END_REF][13] . On the other hand, for a fixed value of λ min , the minimum ∆x verifying the resolution requirements grows roughly linearly with s. As a consequence, for a given value of the maximum frequency to be resolved and for given dimensions of the physical domain of interest, the number of points per direction decreases linearly with s and the total number of nodes diminishes quadratically. This study clearly shows the potential gain in the computational cost offered by the ADESIM sk method over classically-employed FD schemes.

VI. CONCLUSION

An arbitrary-order immersed interface method has been developed for the simulation of the two-dimensional linear propagation of mechanical waves in heterogeneous media including fluid and solid domains. The proposed technique uses ADER schemes of arbitrary accuracy in space and in time to solve the governing equations and employs the "explicit simplified interface method" (ESIM) by Lombard and Piraux [START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF] for the treatment of the interfaces. In order to implement the ESIM, arbitrary-order spatial derivatives of the interface conditions are required. To this end, an algorithm for their numerical computation, which does not require the knowledge of their explicit analytical expressions, has been derived. This procedure could be extremely useful for future developments in three dimensions. Two numerical experiments have been finally examined: the propagation of a plane wave through a flat fluid/solid interface and the propagation of a cylindrical pulse through an aluminum shell immersed in water and filled with water. The improvement of the convergence rate and of the numerical accuracy when increasing both the order of the ADER scheme and of the interface treatment has been more specifically demonstrated. In particular, a sixthorder ADER method coupled with a fifth-order ESIM allows achieving a resolution of about 7.6 points-per-wavelength, which is comparable to the accuracy of very-high-quality solvers, such as the well-known SPECFEM code.

The present study leaves various unanswered questions that will be the subject of future investigations. Firstly, specific strategies must be developed for the treatment of non-smooth and crossing interfaces. Secondly, a stability analysis is required to rigorously define the range of applicability of the present numerical technique. Finally, the method must be extended to three-dimensional configurations in order to investigate more realistically the propagation of mechanical waves in heterogeneous media.

FIG. 4 . 4 FIG. 5 . 2 - 2 , 6 .

 445226 FIG. 4. (color online) First test case. Stress fields σ yy at three different instants of time obtained through the scheme ADESIM 6 -5.

FIG. 6 .

 6 FIG. 6. (color online) First test case. (a) Reflected pressures at the point (-1, 0) m. (b) Transmitted P-waves and (c) transmitted S-waves at the point (0.8, 0) m. (d)-(e)-(f) Corresponding one-sided energy spectral densities. Results obtained analytically (black lines) and numerically, with ∆x = 1/128 m, through the following methods: ADESIM 2 -2 (purple lines), ADESIM 4 -3 (blue lines), ADESIM 6 -3 (green lines), ADESIM 6 -5 (orange lines).

  FIG. 7. (color online) Second test case. Stress fields σ xx at four different instants of time obtained through the scheme ADESIM 6 -5.

FIG. 8 .

 8 FIG. 8. (color online) Second test case. (a)-(b) Signals at the point of coordinates (-0.5, 0) m obtained analytically (black lines) and through the schemes ADESIM 2 -2 (purple lines), ADESIM 4-3 (blue lines), ADESIM 6-3 (green lines), ADESIM 6-5 (orange lines). (c) Corresponding onesided energy spectral densities. (d) Relative errors for the absolute values of the Fourier transforms of the numerical signals.

TABLE I .

 I Second test case. Numerical accuracies in terms of points-per-wavelength.

	ADESIM 2 -2 ADESIM 4 -3 ADESIM 6 -3 ADESIM 6 -5	15.38 34.39 51.75 52.63	26.01 11.63 07.73 07.60

Method sk Max. resolved f (Hz) λ min /∆x

ACKNOWLEDGMENTS

This work was performed under the Project SAC000D04-High-Resolution Low-Frequency Synthetic Aperture Sonar (HRLFSAS) of the STO-CMRE Programme of Work, funded by the NATO Allied Command Transformation.