Giuseppe Rega
email: giuseppe.rega@uniroma1.it

Stefano Lenci
email: lenci@univpm.it

Nonlinear Mechanical Oscillators

Keywords: Dynamical integrity measure, safe basin erosion, nonlinear mechanical oscillators, control, reduction of erosion, escape

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

This work overviews and continues recent investigations by the authors [START_REF] Rega | Identifying, evaluating, and controlling dynamical integrity measures in nonlinear mechanical oscillators[END_REF] on the dynamical integrity of nonlinear mechanical oscillators, an issue which was formerly addressed by Thompson and co-workers (Thompson, 19891 Soliman and Thompson, 19891 Lansbury et al., 1992).

It is nowadays realized that, for safe practical applications, dynamic attractors must be paralleled by uncorrupted basins. Indeed, eroded basins constitute a critical state corresponding to impending system failure, which can be instantaneous and dramatic in the presence of perturbations or uncertainties. Yet, dealing with a system only when the erosion is totally prevented may be too conservative, because it can actually survive safely well above the relevant threshold, if the erosion is not sharp.

These considerations call for a detailed investigation of the issue of dynamical integrity, which is the first objective of this work and schematically consists of various basic steps.

1. Choosing the right definition of "safe basin," i.e., of what is dynamically acceptable. 2. Choosing an appropriate integrity measure quantifying the magnitude of the safe basin.

Both points strongly depend on the problem at hand and may significantly vary from case to case. 3. Investigating the evolution of the basin for a varying problem parameter. When the integrity is reduced throughout this process, we name it "erosion". From an engineering viewpoint this phenomenon can be studied by the erosion profiles, which provide key information such as, the start and end points of erosion. As a matter of fact, they also furnish fundamental insight to the design of the controllable parameter, allowing for its optimal selection if the erosion profile has maxima, or suggesting to stay far away from relevant dangerous sharp falls. 4. The last point consists of analyzing the outcomes of safe basin collapses, which correspond to the onset of different, system-dependent, phenomena, which may or may not represent the actual collapse of the structure, and can therefore have different degrees of dangerousness.

All of these analyses should be paralleled by theoretical work aimed at studying the topological mechanisms responsible for integrity reduction and for the bifurcation event determining the final destruction of the safe basin. However, while providing fundamental understanding, these investigations are not likely to provide meaningful practical information, so they will be deferred to future investigations.

The previously mentioned topics concern the analysis of given systems and the understanding of their dynamical behaviour. Of course, from an engineering point of view it can be useful, or even necessary, to act in such a way to maintain adequate levels of dynamical integrity within fairly large parameters regions. Thus, control techniques able to provide the required performance are welcome and worthy of application, and permit better practical use of the oscillator.

In fact, the second objective of this work is to determine whether it is possible to preserve the dynamical integrity under severe working conditions, or, equivalently, to reduce/postpone the safe basin erosion. A technique previously developed by the authors (Lenci and Rega, 2004a) for controlling nonlinear dynamics and chaos is applied. The method consists of optimally eliminating appropriate homo/heteroclinic bifurcations governing system dynamics by varying the shape of the excitation, and relies on the observation that the erosion is usually triggered by the global bifurcation of a given (usually hilltop) saddle, although it can also involve other topological events.

Different systems (e.g., smooth versus non-smooth, hardening versus softening, symmetric versus asymmetric, etc.), chosen to cover some main mechanical, dynamical and control aspects, are addressed and compared with each other as regards the integrity behavior, the features of erosion and final escape, and the performances of the control procedure. Both classical (e.g., Helmholtz, Duffing, the rocking block) and non-classical (a model of an electrodynamically actuated microbeam) oscillators, representing actually discrete systems or unimodal approximations of continuous structures, are considered. Uncontrolled and controlled erosion profiles are built, and investigated in detail.

While the various mechanical systems have been separately, and independently, investigated in detail in previous authors' works (Lenci and Rega, 2003a1 2003b1 2005a1 2006), they are herein addressed and overviewed in a unitary framework. As a matter of fact, apart from some minor specific points, the original part of the work is a detailed and systematic comparison amongst the various systems, aimed (i) at enlightening on how the dynamical integrity is influenced by their mechanical properties, in general, and (ii) at illustrating the effects of control in terms of dynamical integrity, in particular. The latter aspect is particularly worthy, because it permits to compare the practical effects of control, and to judge on its effectiveness with respect to a large number of mechanical, bifurcational and dynamical issues, thus complementing and completing on the technical side the more theoretical comparison performed elsewhere (Lenci and Rega, 2004a). The joint results of these two works provide the necessary background and the overall knowledge necessary to acquire reliability on control.

The paper is organized as follows. In Section 2 the definitions of safe basin, the evaluation of integrity through appropriate measures, the reduction of erosion by varying systems parameters, and the maintenance of integrity by control are briefly addressed. The four reference mechanical systems are studied in Section 3, while Section 4 is entirely devoted to a detailed comparison and discussion of the principal properties of dynamical integrity (including the onset of out-of-well attractors after escape) and control, as highlighted by the considered sample systems. The paper ends with some conclusions and suggestions for further developments (Section 5).

DYNAMICAL INTEGRITY: DEFINITIONS, EVALUATION, REDUCTION (EROSION) AND CONTROL

The fundamental issues necessary for an accurate investigation of the dynamical integrity, of its evaluation, of its erosion and of its control have been studied in [START_REF] Rega | Identifying, evaluating, and controlling dynamical integrity measures in nonlinear mechanical oscillators[END_REF], a paper which the reader is referred to for details. The main aspects of that work are summarized in the following subsections to set the proper research framework.

Safe basin

The safe basin is the union, in phase space, of all initial conditions guaranteeing some specific, dynamically acceptable, response performances, which are usually the convergence in time toward one or more attractors, or the non-escape from a potential well. It can be the safe basin of a specific attractor, in this case coinciding with a classical basin of attraction [START_REF] Hsu | Cell-to-Cell Mapping: A Method of Global Analysis for Nonlinear Systems[END_REF], or it can refer to a potential well, in this case being the union of the basins of all attractors belonging to the considered well.

The safe basin can be built either by considering only the steady state dynamics, which makes sense when one is interested in the system stationary regime, or by paying attention also to the transient dynamics, in cases where the short-term dynamics is the most important part of the response or some unsafe phenomena -such as a temporary escape from the potential well -may occur during the transient. Technically, they differ because the former basins are bounded by the stable manifolds of some relevant saddles, thus having a clear dynamical meaning, whereas the latter do not possess this property. Yet, in several applications, the ensuing differences are actually minor and this distinction remains mainly a theoretical issue. Accordingly, in this paper only steady dynamics will be considered.

The previous definitions do not care about the phase of the excitation, though it can actually be important in many practical situations. To take this aspect into account, the "true" safe basin is defined as the intersection of the previously defined safe basins when the excitation phase ranges over a time interval of interest (e.g., the period for periodic excitations, or the period of free vibration for impulsive excitation). The "true" safe basin is the smallest phase-independent set of initial conditions leading to "safe" dynamics.

Integrity measures

Various parameters have been proposed to measure the magnitude of the safe basins (Soliman and Thompson 19891 Rega and Lenci, 2005). They have different meaning, relevance and properties. Only two of them will be used in this work.

The Global Integrity Measure (GIM) is the normalized hyper-volume (area in 2D cases) of the safe basin. It is the most intuitive and easy integrity measure, but it is not satisfactory in all cases in which the safe basin is strongly fractal. Normalized means that the hyper-volume of the safe basin corresponding to the actual value of a varying parameter is divided by that corresponding to a reference value, so that GIM is a dimensionless number.

The Integrity Factor (IF) is the normalized radius of the largest hyper-sphere (circle in 2D cases) entirely belonging to the safe basin. The IF is as computationally easy as the GIM, but it succeeds in eliminating the unsafe fractal tongues from the integrity evaluation. In fact, it is a measure of the compact part of the safe basin, which is the solely region which guarantees the dynamical integrity of the system.

Reduction of integrity and the erosion profiles

The integrity measures are useful in studying how the system reliability is changed when the parameters vary. Drawing the integrity measure as a function of a varying parameter provides the so-called "erosion profiles", which highlight the variation of dynamic integrity.

In practical applications a very important parameter is the excitation amplitude, and its increase usually determines a reduction of integrity, namely, a loss of structural safety, which justifies the name erosion. The erosion profiles for various mechanical oscillators will be reported in Section 3.

The erosion can be studied through the detection of proper invariant manifolds, which are the basin boundaries, and the analysis of their evolution, in particular by studying the relevant global bifurcations, which are the key mechanisms for erosion. This is exact if the safe basins coincide with classical basins of attraction, whereas only approximate, but in any case enlightening, in the other cases.

The reduction of integrity is commonly triggered by homoclinic or heteroclinic bifurcations allowing the basins of surrounding attractors to enter the safe basin and to reduce their magnitude. Then, it proceeds with different topological mechanisms, which are extremely complex and may involve secondary homo/heteroclinic bifurcations, as already mentioned. The erosion ends with the complete destruction of the safe basin and with the onset of other dynamical regimes, which may be desirable (e.g., in the case of anti-control, in MEM switches), simply tolerable (e.g., cross-well dynamics in hardening systems) or destructive (ship capsize, rigid block overturning, MEM sensors).

Erosion reduction by control

The basin erosion clearly reduces the dynamical integrity of the system, so that the possibility of its control, namely the increment of the integrity measure, has a great theoretical and practical interest.

With this and other objectives, the authors (Lenci andRega, 2004a, 2004b) have recently developed a general control technique aimed at shifting the homo/heteroclinic bifurcations towards higher excitation amplitudes. Although the idea is general, it has been mainly applied to the bifurcations of homoclinic and/or heteroclinic loops surrounding the potential well(s), those which asymptotically approach the hilltop saddle(s), which can be detected by the Melnikov method (Guckenheimer and Holmes, 19831 Wiggins, 1988).

The global bifurcation can be eliminated in various manners, and in previous works it has been done by optimally varying the excitation shape, i.e., by adding controlling superharmonics to a given harmonic term, to be considered as the reference excitation.

Two different approaches have been proposed. If one wishes to shift only a single global bifurcation, then one applies the so-called "one-side" control, which is theoretically very effective (the critical amplitude for bifurcation can grow up to a factor 2) but controls only one part of the phase space. When there are two (or more) relevant global bifurcations, "one-side" control is paralleled by the so-called "global" control aimed at increasing simultaneously all possible homo/heteroclinic bifurcations. While providing a spread control of phase space, this approach is theoretically less performant, and thus it is alternative and not competing with the "one-side" control.

The effectiveness of this control technique in shifting the erosion profiles toward higher excitation amplitudes, which means that the control is effective in reducing erosion, will be illustrated in the following section.

VARIOUS MECHANICAL OSCILLATORS

Four different systems are considered, chosen with the aim of covering the main mechanical, dynamical and control features and investigating how they influence the erosion. A detailed discussion and comparison of various aspects of the dynamics exhibited by the four oscillators will be the subject of Section 4.

Helmholtz

The dimensionless Helmholtz equation

2 x 3 011 4 x 5 x 3 x 2 6 2 34t5 6 2 1 N 1 j61 2 j 2 1 sin3 j4t 3 6 j 5 (1)
describes a smooth, softening and asymmetric oscillator and governs the dynamics of various softening mechanical systems (Lenci and Rega, 2003a). 2 (4t) is the generic 27/4-periodic external excitation: 2 1 is the amplitude of the basic harmonic, and a measure of the overall excitation amplitude, while 2 j are the amplitudes of the controlling superharmonics, which, together with the phases 6 j , govern the excitation shape and are optimally determined by the control technique (Lenci and Rega, 2003a).

The erosion profiles corresponding to harmonic and "one-side" control ("global" control does not apply because there is a single potential well surrounded by a single homoclinic loop) with either one or two controlling superharmonics are reported in Figure 1. Here the safe basins are the classical basins of (in-well) bounded attractors, and the integrity is measured by the GIM.

The excitation frequencies considered are close to the vertex of the escape triangle in the excitation parameters plane [START_REF] Szemplinska-Stupnicka | Cross-well chaos and escape phenomena in driven oscillators[END_REF]. Thus, we are considering the most dangerous frequencies (those for which escape occurs for the lowest value of excitation amplitude), where even small improvements are very useful.

The vertical segments correspond to the numerically detected homoclinic bifurcation of the hilltop saddle, which, as expected, is clearly seen to be the dynamic phenomenon triggering the erosion in all the cases considered.

Figure 1 clearly shows how the control is effective in reducing erosion. An in-depth discussion of Figure 1 is, however, postponed to Section 4, where the diagrams will also be compared with those of other mechanical systems.

Duffing

The dimensionless Duffing equation

2 x 3 011 4 x 5 x 2 3 x 3 2 6 2 34t5 (2)
describes a smooth, hardening and symmetric oscillator and governs the dynamics of various hardening mechanical systems (Lenci and Rega, 2003b). The excitation 2 (4t) is the same as in (1).

The erosion profiles corresponding to harmonic, "one-side" and "global" (which now applies) controls with a single controlling superharmonic are reported in Figure 2 for 4 6 018, a value close to the vertex of the V-shaped region of cross-well chaos [START_REF] Szemplinska-Stupnicka | Cross-well chaos and escape phenomena in driven oscillators[END_REF] and therefore practically important. Here the safe basins are still the classical basins of attraction of in-well attractors, but the integrity is now measured by the IF. Figure 2 shows how the global optimal excitation (curve (b)) is able to shift the erosion curve toward larger amplitudes with respect to the harmonic excitation (curve (a)), namely, it is effective in reducing erosion. The (right) one-side control has different properties. In fact, due to the asymmetry of this excitation (Lenci and Rega, 2003b), the dynamics in the left and right wells are different: the erosion curve (c) l of the left uncontrolled well is lower than curve (a)1 namely, there is strong erosion in the uncontrolled potential well. In the right controlled well, on the other hand, basically there is no erosion at all (curve (c) r ), thus showing the local character of the one-side control.

Both one-side erosion curves end at 2 1 7 01066, where the last attractor belonging to the left potential well disappears by a saddle-node bifurcation (Lenci and Rega, 2003b). For 2 1 8 01066 (and up to 2 1 7 010951) there are only confined attractors in the right controlled well, so that their safe basin coincides with the whole phase space and erosion no longer makes sense. We close this section by noting that the curve (c) r does not fall down, in agreement with the fact that the right homoclinic bifurcation occurs at 2 r 19cr 6 010807, a value well above the point where the curve disappears.

Rigid block

The nonlinear dynamics of a slender rigid block rocking around base corners in a constant gravitational field is governed by the well-known Housner model (Lenci andRega, 2005a, 2005b):

2 x 3 4
x 5 x 3 012 3 2 34t5 6 09 x 8 08 2 x 3 4

x 5 x 5 012 3 2 34t5 6 09 x 08

4 x3t 3 5 6 r 4 x3t 5 59 x 6 01 (3)
In ( 3) is the damping coefficient, r the restitution coefficient, and the excitation, representing the dimensionless horizontal acceleration of the foundation, is the same as in (1) with 6 j 6 6 j 3 72.

Contrary to previous examples, for the rigid block there is no natural frequency around which to focus numerical analyses. Thus, there are no preferred values of the excitation frequency, and we choose 4 6 315. For this value and for 6 0102 and r 6 0195, the four erosion profiles with harmonic excitation and corresponding to classical and "true" in-well safe basins as well as to GIM and IF are reported in Figure 3, and are discussed in Section 4. Note that the vertical line corresponds to the homoclinic bifurcation threshold 2 19cr 6 01094 and that the curves are normalized with respect to their values at 2 1 6 0105.

Figure 3 aims at comparing two safe basin definitions and two integrity measures. To visualize the differences between classical and "true" in-well safe basins we report two relevant examples in Figure 4. Not only are the magnitudes strongly different, but the degree of fractality is also much reduced in the "true" safe basin, which is more compact and closer to the interior circle involved in the definition of IF, which is also reported.

In Figure 3 we fixed the kind of excitation and varied both the safe basin definition and the integrity measure. The opposite is done in Figure 5, where we report the erosion profiles for the classical safe basin of the potential well and the IF, for both harmonic and "global" control with a single added superharmonic excitation ("one-side" control is herein theoretically expected not to work (Lenci and Rega, 2005a)). Once more, control is shown to be effective in shifting the erosion profile toward higher excitations, in particular the sudden fall occurring just ahead of the homoclinic bifurcation, although it becomes ineffective afterwards.

MEMS

The dimensionless equation 2 x 3 0101 4

x 3 x 3 x 3 5 0114 3 2 34t5 31 5 x5 2 6 0 ( 4 ) describes the single d.o.f. nonlinear dynamics of a thermoelastic rectangular microbeam subjected to an axial load and to a concentrated electrodynamic transverse force applied at midspan, fixed at both ends and supposed to be in an ultra-high vacuum environment [START_REF] Lenci | Control of pull-in dynamics in a nonlinear thermoelastic electrically actuated microbeam[END_REF]. 2 (4t) is still given in (1), and here it is proportional to the applied oscillating voltage.

The effects of control on the erosion profiles are shown in Figure 6 for both IF and GIM. Here the classical safe basin of the potential well is used, and only one "one-side" controlling superharmonic is added to the basic harmonic excitation ("global" control does not apply). As for the Helmholtz and Duffing oscillators, the frequency 4 6 017 is close to the vertex of the V-shaped region of escape [START_REF] Lenci | Control of pull-in dynamics in a nonlinear thermoelastic electrically actuated microbeam[END_REF].

The three curves in Figure 6 correspond to harmonic, generic control and optimal control excitations. It is evident that the control is able to shift the erosion profiles, and that the shift is larger with the theoretically determined optimal value of the amplitude of the controlling superharmonic.

To better highlight the effect of the increasing amplitude of the controlling superharmonic, we report in Figure 7 the erosion profiles for the fixed excitation amplitude 2 1 6 010025, which corresponds to the sudden fall of the curve of harmonic excitation in Figure 6, Figure 6. MEMS: Erosion profiles for harmonic (2 2 2 1 6 0), control (2 2 2 1 6 015) and optimal control (2 2 2 1 6 116591) excitations for 4 6 017. The homoclinic bifurcation in the case of harmonic excitation occurs at 2 19cr 6 01001078. and for the varying relative amplitude 2 2 /2 1 . Looking at the numerical gains G, which are values of each integrity measure normalized to 1 in correspondence with the harmonic excitation, it is clear that the best result is obtained for 2 2 /2 1 7 116, which is in perfect agreement with the theoretically predicted optimal value 2 2 /2 1 6 116591. This occurs with both IF and GIM, although the former integrity measure is more effective in highlighting the benefits of control. The possibility to satisfactorily also apply non-optimal control is discussed at the end of Section 4.

DISCUSSION OF UNCONTROLLED/CONTROLLED DYNAMICS OF DIFFERENT OSCILLATORS

This section is devoted to a discussion and comparison of dynamical integrity features of various mechanical oscillators, as briefly illustrated in Section 3. The basic observation is that control is always able to shift the erosion profiles corresponding to reference harmonic excitation towards higher excitation amplitudes (Figures 1, 2, 5 and 6), and thus it is broadly effective, although to a different extent, in increasing the range of dynamic reliability of the oscillator, this representing a practically appealing feature.

The other points of specific comparison are listed below by grouping them within two broad classes of mechanical/bifurcational issues and integrity/control issues. Prior to proceeding further, it is worth noting that due to the large variety of mechanical systems considered, the properties reported below are suggested to be very general, but of course they hold rigorously only for the cases considered, and different situations can in principle occur for other values of the parameters and/or for different oscillators.

Mechanical and bifurcational issues

Hardening (Duffing) versus softening (Helmholtz, rigid block, MEMS) systems. The main difference between hardening and softening systems is related to the out-of-well attractor onto which the system settles down at the end of the erosion process. It is a cross-well attractor for the Duffing oscillator, which in practical applications may be undesirable but not destructive, whereas it is an attractor corresponding to the actual system failure for all of the others.

In terms of erosion, the difference is that for (multi-well) hardening systems the eroding attractor(s) is (are) that (those) of the neighboring potential well, and vice versa, and we observe mutual tangling of the relevant safe basins. For softening systems, on the other hand, the safe basin of the single well is always eroded by the out-of-well, destroying, attractor. This difference, however, stays in the background, and does not emerge in the erosion profiles. In fact the sole difference occurs in the erosion curves of the Duffing oscillator with "one-side" control, which are shown to be interrupted by the disappearance of the eroding attractor, a fact that is likely non-generic and that cannot occur for softening systems where the destroying attractor always exists.

As there are no more specific differences attributable to the different mechanical nature of the systems, we can conclude that the control effectiveness is independent of the kind of nonlinear behavior.

Symmetric (Duffing, rigid block) versus asymmetric (Helmholtz, MEMS) systems. Another mechanical characteristic that is illustrated by the four systems considered is the symmetry of the restoring force. The effects of the symmetry are that (i) in the Duffing oscillator the two homoclinic bifurcations occur simultaneously for both harmonic and "global" control excitations (but not for the "one-side" control, which is an asymmetric excitation making the system dynamical behavior asymmetric), and that (ii) in the rigid block the potential well is surrounded by a couple of heteroclinic orbits, which become a unique, structurally stable, homoclinic loop in the presence of imperfections, thus rendering the "one-side" control inef-fective in this case (see the following point "Homoclinic versus heteroclinic bifurcations"). While these properties are quite important for the theoretical developments of control, they appear to have minor consequences on the dynamical integrity.

In fact, the results of Section 3 show that, in terms of erosion profiles, there are no significant differences between symmetric and asymmetric systems, with both harmonic and control excitations, so that this mechanical feature does not seem to influence the dynamical integrity behavior.

Among the others, this conclusion is likely to be particularly related to the specific cases considered in Section 3. In fact, it can be conjectured that in other asymmetric systems, such as, for example, the Helmholtz-Duffing oscillator, the asymmetry may play a non-negligible role in terms of dynamical integrity, as suggested by the theoretical developments in Lenci and Rega (2004b).

Smooth (Helmholtz, Duffing, MEMS) versus non-smooth (rigid block) systems. The third mechanical property that can be discussed on the basis of the sample oscillators of Section 3 is the system smoothness versus non-smoothness. This difference has meaningful consequences in the analytical detection of the homo/heteroclinic bifurcations, which is, however, a technical point, while it does not seem to have consequences in terms of dynamical integrity and erosion profiles. In fact, even the perfect coincidence between the global bifurcation and the sudden fall of the profiles (see the following specific point), which is seen to occur only in Figure 5, cannot be attributed to the system non-smoothness, as it does not occur, for example, in Figure 3.

More generally, based on our experience, we can say that the piecewise nature of the system does not seem to entail specific characteristics in terms of either dynamical integrity or of its control. This conclusion is further supported by the results of previous numerical investigations of the nonlinear dynamics and control of an inverted pendulum with rigid lateral barriers [START_REF] Lenci | Controlling nonlinear dynamics in a two-well impact system. Parts I and II[END_REF], herein not reported, where no specific outcomes that could be clearly ascribed to the presence of a piecewise restoring force were highlighted. Of course, this does not mean that any specific effects may actually exist as regards the bifurcational behavior, but just that they are not robust enough to systematically affect the integrity scenario. We conjecture that there might be pathological effects requiring a detailed investigation.

Global bifurcations triggering the erosion.

A characteristic observed in all cases apart from the MEMS (Figure 6) is that the starting points of erosion are just ahead, and in good agreement with, the homo/heteroclinic bifurcations of the hilltop saddles. Moreover, in the case of the rigid block (Figure 5), this agreement is indeed a perfect coincidence.

This property confirms that these global bifurcations are the necessary pre-requisite for erosion, which successively evolves through other, topologically more complex, bifurcational events [START_REF] Thompson | Chaotic phenomena triggering the escape from a potential well[END_REF]) up to the system escape. Along with the associated possibility of analytical treatments, this justifies why the control method, which in principle can be applied to any global bifurcation (Lenci and Rega, 2003c), has been correctly applied to the homo/heteroclinic bifurcations of the hilltop saddles.

In the case of MEMS, the control is in any case effective, in spite of the fact that the sudden fall in the profiles is likely due to a secondary heteroclinic bifurcation embedded in the dynamics and not to the homoclinic bifurcation of the hilltop saddle [START_REF] Lenci | Control of pull-in dynamics in a nonlinear thermoelastic electrically actuated microbeam[END_REF]. This further underlines the robustness of control, which shows to be practically effective also beyond theoretical expectations. Homoclinic (Helmholtz, Duffing, MEMS) versus heteroclinic (rigid block) bifurcations. The cases of Section 3 permit one to compare systems with a couple of heteroclinic orbits surrounding the potential well (rigid block) with systems with homoclinic orbits surrounding the unique (Helmholtz, MEMS) or the two (Duffing) potential wells. This difference has consequences in terms of theoretical implementation of control, and it is at the base of the non-effectiveness of "one-side" control for the rigid block. In fact, if, by "one-side" control, we eliminate only one heteroclinic bifurcation, say, the upper, we inhibit the penetration of fractal tongues from the upper part of the phase space, but we allow for their penetration from the lower part, where the other heteroclinic bifurcation not only has not been removed, but has been amplified (Lenci andRega, 2005a, 2005b).

Apart from this aspect, which has also consequences in terms of maximum expected theoretical performance (which is larger with the "one-side" control (Lenci and Rega, 2004a)), there are no other apparent consequences of this dynamical feature in terms of erosion, as, for example, the erosion profiles of Figure 5 are similar to those of the other mechanical oscillators.

Integrity and control issues

Practical and theoretical saved regions. Figures 1, 2, 5 and 6 show that there is a welldefined interval of amplitudes, ahead of the fall in the profiles with harmonic excitations, where the controlled profiles are above the harmonic profile and thus the control is effective in maintaining the dynamical integrity. This zone, which has been named the practical saved (or controlled) region in previous authors' works, permits one to appreciate, even visually, the advantages of using control excitations.

The practical saved region is usually comparable with the theoretical saved region, which is in between the homoclinic bifurcation thresholds for harmonic and control excitations, and this proves how theoretical results are indeed practically realized.

Furthermore, according to the analytical predictions, which show that the theoretical saved region is "large" (increment of the critical threshold of 41% with one superharmonic) for "one-side" control and "small" (increment of 15%) for "global" control (Lenci and Rega, 2004a), it is observed that also the practical saved region is larger in the case of "one-side" control (Figure 1) and smaller -but in any case present -with "global" control (Figures 2 and5).

Various "failure" phenomena for softening oscillators. In softening systems, the onset of the out-of-well attractor is usually destructive, and leads to structural failure. This corresponds to different physical situations: the ship capsizing for the Helmholtz oscillator, overturning for the rigid block, dynamic pull-in for MEM sensors. However, it is worth noting that there exist softening mechanical oscillators for which the onset of out-of-well dynamics is not destructive and, on the contrary, is desirable. An example is given by the pull-in solutions in MEM switches [START_REF] Nayfeh | Dynamic pull-in phenomenon in MEMS resonators[END_REF].

While the shift of the erosion toward higher excitation amplitudes is highly positive because it increases the range of system dynamic reliability, the onset of the failure attractor is not strictly related to the erosion behavior, but rather to the post-erosion features: the principal two are discussed in the next few points.

Delay of escape.

While being effective in shifting the erosion profiles, the control is not expected to be able to also shift the escape, i.e., the end of the erosion profiles. This surprisingly occurs for the Helmholtz oscillator in the range 4 6 0181 5 0185 (Figures 1(b) and 1(c)), i.e., just across the vertex of the V-shaped region of escape, where even the slight increment of escape threshold in Figures 1(b) and 1(c) is very useful (Lenci and Rega, 2003a). This highlights an unexpected practical performance of control, which hopefully holds also around the vertex of the V-shaped regions of the other softening oscillators, although this has not been directly checked.

In the other cases reported in Section 3, on the other hand, the control does not affect the escape threshold (Figures 1(a), 1(d) and6) or even reduces it, as it is evident in Figure 5, and also in Figure 2. However, it is worth noting that this is not actually a drawback of the control method, whose implementation has not been designed to increase the escape threshold, which could be done in principle and would lead to different optimal excitations, but just to shift the erosion profiles -rigorously, to increase the triggering point of erosion represented by the homo/heteroclinic bifurcation of the hilltop saddle. And the control has shown to be fully effective with respect to this goal.

Sharpness of the profiles.

A sudden fall (the so-called "Dover cliff" [START_REF] Lansbury | Basin erosion in the twin-well Duffing oscillator: Two distinct bifurcation scenarios[END_REF]) is observed in most erosion profiles apart from those of the rigid block (Figure 3, which shows only a relatively smooth fall, and Figure 5, where the sudden fall is marked with the control yet less marked with harmonic excitation), and from that of Figure 6(b) (which indeed would also show a -more marked -fall if drawn on a different scale). This shows that the out-of-well eroding tongues rapidly penetrate and erode the safe basin. An explanation of the sharpness in terms of the influence of resonance on the erosion process is reported next.

Sharpness is unpleasant from a practical point of view, because the sudden loss of integrity means that the system rapidly becomes unsafe, and small imperfections or uncertainties can lead to instantaneous failure of the structure. This is further supported by the fact that the basin boundaries become fractal when approaching the peak of the fall.

While shifting the profiles forward, the control does not show to be able to reduce sharpness, which instead is slightly increased in some cases (e.g., Figure . 1(a) and 1(b)). Thus, while enlarging the range of system reliability (by the saved region), the control is not able to eliminate its final sensitivity to small perturbations, which can be interpreted as a sort of practical structural instability.

Erosion of systems without (rigid block) and with (Helmholtz, Duffing, MEMS) internal frequency. When compared to the dullness of the diagrams in Figure 3, the sharpness of the other profiles can be attributed to the influence of the resonance on the erosion process, or better, on the V-shaped region of escape [START_REF] Szemplinska-Stupnicka | Cross-well chaos and escape phenomena in driven oscillators[END_REF], which is a direct consequence of a resonant excitation. In fact, the rigid block is the only system considered with no internal frequency, whose effect, whether existing and externally excited, is of course to enhance the high amplitude response features ending in escape. This conjecture is further supported by the observation that also in the Helmholtz oscillator, far enough from the vertex of the escape region (e.g., Figure 1(a)) and therefore far away from resonance, the profiles are less sharp than those close to the vertex (e.g., Figure 1(b) and 1(c)).

More specifically, we note that in the neighborhood of the vertex of the V-shaped region the critical threshold for escape, where the erosion ends, is strongly reduced (Szemplinska-Stupnicka, 19921 Lenci and Rega, 2003a) and approaches the homo/heteroclinic bifurcation threshold, where the erosion starts, so that erosion must necessarily develop suddenly. When this closeness is lost, either due to lack of resonance or far away from the resonant frequency, the erosion can, and actually does, develop gradually.

" One-side" versus "global" control. The comparison between "one-side" and "global" controls on a single system can be done only for the Duffing oscillator (Figure 2), which is the unique system with two potential wells. In fact, in the Helmholtz oscillator and in the MEMS there is just one homoclinic bifurcation, so that "global" control does not apply. In turn, for the rigid block, there is a single potential well surrounded by a couple of heteroclinic connections, so that "one-side" control, which eliminates only one of the two bifurcations and amplifies the other (Lenci andRega, 2005a, 2005b), cannot be effective (as already noticed in the point "Homoclinic versus heteroclinic bifurcations") and indeed it has not been used.

Unfortunately, the comparison in Figure 2 is not fully enlightening because in the case considered the erosion curves of "one-side" control end well below the relevant homoclinic bifurcation, so that the extent of the corresponding practical saved region cannot be numerically evidenced. We can only say that (i) the erosion profile of "one-side" control, whether existing, is certainly above the other profiles, according to the theoretical predictions, which suggest the "one-side" control performs better than the "global" control and, of course, harmonic excitation, and that (ii) on the basis of Figure 2 it is possible to conjecture, by looking at the final part of curve (c) r , that if it were hypothetically continued, it would remain above curve (b) up to a certain limit.

The comparison between "global" control of the Duffing system and "one-side" control of the other systems, which is of course only slightly meaningful, seems to confirm the better practical performances of "one-side" control (in this respect, see also the point "Practical and theoretical saved regions").

Classical versus true safe basins and GIM vs IF with harmonic excitation. This question is basically addressed in Figure 3. The most important conclusion which can be drawn on the basis of this picture is that classical safe basins underestimate erosion, which is instead better captured by "true" safe basins. The differences are dramatic when the excitation amplitude increases, as also confirmed by Figure 4.

Still referring to classical safe basins, we see that IF and GIM basically coincide up to 2 1 7 0115, where they diverge in correspondence with a step likely due to a secondary global bifurcation. This shows how GIM performs less well than IF in detecting system evolutions.

In the case of "true" safe basins, an opposite situation is observed. Here, IF is systematically greater than GIM, the differences being not negligible. Thus, one can argue that in this case the GIM is better than the IF, but this requires more confirmation and we feel that it is not a general property. However, both measures correctly capture the global phenomena and provide qualitatively similar erosion profiles.

The comparison between the two IFs and the two GIMs shows that while the former are relatively close to each other, the latter are very far apart. Thus, the IF (GIM) seems to be a measure low-(high-) dependent on the excitation phase. This property shows that (i) the former measure is intrinsically better than the latter, and that (ii) if use is made of IF instead of GIM, we automatically "rule out" -or better strongly reduce -the effects of the excitation phase without using the "true" safe basins, whose construction is computationally more expensive, with a great operational advantage.

GIM versus IF for harmonic and control excitation. A detailed comparison between IF and GIM as obtained with harmonic, (non optimal) control, and optimal control excitations is reported in Figure 6. First, we note, again, that IF and GIM basically coincide up to the vertical step due to a secondary bifurcation. The overall curves scenario (i) confirms the findings for the rigid block under harmonic excitation (see previous point) and extends them to MEMS, (ii) also extends them to control excitations, and (iii) further supports the conclusion that IF is better that GIM in detecting system evolutions, as confirmed by the fact that IF highlights the erosion step better than GIM and that, afterwards, IF is always smaller than GIM for all excitations.

Both IF and GIM show the improvements obtained on passing from harmonic to (nonoptimal) control excitation, and also highlight that the largest shift is obtained with optimal control excitation. More specifically, the horizontal shifts of the profiles are comparable, while the vertical shift for a fixed excitation amplitude is more marked with IF, due to sharpness of the relevant profiles. This shows once again that IF is a better integrity measure. In any case, the profiles maintain their qualitative shapes on passing from harmonic to control excitations.

Erosion profiles for varying frequency. This point has been investigated for the Helmholtz oscillator (Figure 1) around the critical value 4 6 0181 corresponding to the vertex of the escape triangle [START_REF] Szemplinska-Stupnicka | Cross-well chaos and escape phenomena in driven oscillators[END_REF]. The first property highlighted by Figure 1 is that shift of the erosion profiles is substantially independent of the frequency, showing the robustness of the method and being in agreement with the theoretical predictions, which suggest an 4-independent increment of the critical threshold. The slope of the profiles in the eroding region, on the other hand, varies with the frequency, becoming sharper around the critical frequency for the reasons discussed in "Erosion of system without and with internal frequency".

With control excitations, the shift of the erosion profiles extends just above the relevant homoclinic bifurcation values. In fact, after the beginning of the erosion, the profiles are sharp and rapidly tend to that of the reference harmonic excitation. For "low" frequencies (Figure 1(a)) the three curves cross each other at nearly the same point, after which the GIM of the control excitations is lower than that of the harmonic excitation. Thus, in this frequency range, the control is effective in the saved region, ineffective above it.

Around 4 6 0181 the curves do not cross any more and the escape under control occurs for larger excitation values than with harmonic excitation (Figures 1(b) and (c)1 see also the point "Delay of escape"). Here the initial sharpness of the harmonic profile is slightly increased by control, thus amplifying the sort of instability discussed at the point "Sharpness of the profiles".

For "large" frequencies (Figure 1(d)), the profile in the case of control with one superharmonic tends to recover the harmonic profile, while with two superharmonics it rapidly tends to zero, and the escape threshold decreases. Apart from this point, Figure 1(d) seems to confirm that also for large frequencies the control is effective in the saved region, as theoretically expected, while it is substantially ineffective above it.

Although control is not expected to improve the escape behavior, it is worth noting that for low and high frequencies the control escape threshold is decreased with respect to the harmonic one, while the opposite occurs around the critical frequency. Thus, the V-shaped escape region is both opened and raised (Lenci and Rega, 2003a). The fact that for small and large frequencies the escape is lowered by control is not actually dangerous, because in the lateral ranges the escape already occurs for fairly large amplitudes with harmonic excitation, so that control is practically unnecessary.

Effects of increasing number of controlling superharmonics. The effect of the increasing number N of controlling superharmonics has been investigated in the case of Helmholtz oscillator (Figure 1), and therefore with only "one-side" control.

The general conclusion is that the larger the number of superharmonics, the larger the horizontal shift of the profiles (apart from the pathological case of Figure 1(d)) as well as the practical saved region, in agreement with the theoretical predictions. In Figure 1(a) this correspondence holds up to the end of the saved region, where the control is expected to work, and not further on, where it has no more resources.

A more detailed look at Figure 1 shows that the relative increment due to the addition of the second superharmonic is lesser than the increment due to the first superharmonic. This is also in qualitative agreement with the theoretical predictions (Lenci and Rega, 2003a), which suggest the increment from harmonic to control with N 6 1 to be about 1.41, while the increment from harmonic to control with N 6 2 is about 1.62, so that the first one is expected to be about 40% and the successive relative increment is expected to be 11621141 6 11148, namely, about one third of the first increment (compare with Figure 1(a)).

Effects of relative amplitudes of added superharmonics. This point, which indeed entails checking the practical effects of non optimal control superharmonic terms, has been investigated with reference to MEMS.

The effects of the increasing excitation amplitude 2 1 with a non optimal fixed value of the superharmonic relative amplitude 2 2 /2 1 are shown in Figure 6 and have been illustrated in other points. The opposite case of fixed excitation amplitude and varying 2 2 /2 1 is illustrated in Figure 7, which reports IF and GIM divided by the respective values of integrity measure with the harmonic excitation. Thus, the curves in Figure 7 measure the practical increments with respect to the reference harmonic excitation.

The most important point is that the maximum of both curves is attained for the same value of 2 2 /2 1 , which actually coincides with the theoretically predicted value [START_REF] Lenci | Control of pull-in dynamics in a nonlinear thermoelastic electrically actuated microbeam[END_REF]. We can then conclude that theoretical predictions are well verified in practice, somewhat beyond expectations, because the controlled homoclinic bifurcation only triggers the erosion and is not directly involved in the subsequent basin evolution, so that some major discrepancies would not actually be surprising.

While sharing the same qualitative shapes, the two curves of Figure 7 have important quantitative differences due to their different ability to measure the effect of penetrating eroding tongues on the integrity, which is also the reason why IF is usually sharper than GIM (see the appropriate point). In particular, with IF, the optimal increment with respect to the harmonic excitation is about 2.22, which is significantly larger than the theoretical value 1.41 (see previous point), thus showing how practical performances are somehow better than theoretical expectations1 instead, with GIM, the optimal increment is significantly smaller, being about 1.08. Yet, it is still larger than 1, thus showing how also this measure is capable of highlighting the benefits of control, though to a reduced extent. These results suggest, again, that the IF is a better measure of integrity in practical applications.

An important conclusion can be drawn from Figures 6 and7, namely that quite good results can be obtained with a superharmonic even smaller than the optimal one, at least for the case considered. This is due to the fact that, for example, for 2 2 2 1 6 015, the theoretical [START_REF] Lenci | Control of pull-in dynamics in a nonlinear thermoelastic electrically actuated microbeam[END_REF] and practical gains are still larger than 1, though not optimal. This may be useful in applications, because the optimal excitation may be too expensive, and one may prefer cheaper control excitations providing still satisfactory, although not optimal, results.

Control effort. An issue which cannot be judged based on the results summarized in Section 3 but is important is the control effort needed to achieve the enhancement of the integrity. To quantify this aspect, one can naturally refer to the relative magnitudes 2 i /2 1 of the added controlling superharmonics, which also have the advantage of being dimensionless numbers. Other choices are of course possible, and the smaller 2 i /2 1 the cheaper the control, even if in some cases it might be not so difficult to add high amplitude superharmonics.

The ratios 2 i /2 1 can be expressed in the form 2 i /2 1 6 f 345, where f (.) varies from system to system and is an increasing function of the excitation frequency 4. At first glance this might appear a drawback, which is not actually the case in practice. Indeed, for large values of the excitation frequency, the critical threshold for homo/heteroclinic bifurcation occurs in the various systems at values of the excitation amplitude so high as to need no control. Accordingly, the sub-horizontal plateau of the erosion profiles extends over a quite large range of excitation amplitude values, then falls at the critical thresholds in a less sharp way than in most profiles of Section 3 (see the previous specific points in this section), again confirming the uselessness of control in this case. Control is indeed commonly necessary around resonant frequencies, where the ratios 2 i /2 1 assume reasonable values, i.e., the cost of control is small enough to be proposed in practical applications.

CONCLUSIONS AND FURTHER DEVELOPMENTS

The three steps necessary for correctly evaluating the dynamical integrity of mechanical systems have been critically addressed: definition of the safe basin1 use of an integrity measure characterizing its magnitude1 study of the variation of integrity due to variations of meaningful system parameter(s).

The safe basin has been identified as the union of the basins of attraction of all the attractors belonging to a given potential well, while two different measures, the GIM and the IF have been used. The former is the most natural and simple measure, the latter permits one to rule out the unsafe region of fractality while maintaining a certain algorithmic simplicity.

The variation of integrity when changing the excitation amplitude has been analyzed through the erosion profiles.

The issue of maintaining the dynamical integrity as far as possible by an appropriate control technique has been considered. A control method previously developed by the authors and based on the elimination of the homo/heteroclinic bifurcation of hilltop saddles, which is the event triggering the erosion, has been applied.

Dynamical integrity and control of four different archetypal mechanical systems have been studied for illustrative and comparison purposes, by numerically determining the relevant erosion profiles.

The principal conclusion which can be drawn is that the control is always effective in reducing erosion, to an extent which varies from case to case and which has been investigated in depth. This is proved by the erosion profiles, which always show good agreement with the theoretical predictions. It has been shown that in certain cases the control works well above theoretical expectations.

Several subsidiary specific issues of mechanical, dynamical and control nature have been discussed in detail with the aim of shedding light on the phenomenon of dynamical integrity reduction, and on its control.

As far as further developments are concerned, we mention the study of the erosion of other systems having different mechanical features, such as the Helmholtz-Duffing oscillator, and a detailed study of the erosion of the basin of one single attractor instead of that of the union of in-well attractors. The erosion of out-of-well attractors is also of interest especially in view of specific practical applications concerned with pull-in solutions of MEM switches or rotating solutions of the parametrically excited pendulum.

These are developments of interest from an engineering viewpoint. Also more theoretical developments are worth pursuing. Among others, we quote the in-depth investigation of the topological mechanisms responsible for erosion, in particular as regards the sudden fall in the profiles, along the lines first investigated in [START_REF] Thompson | Chaotic phenomena triggering the escape from a potential well[END_REF].

Finally, it would certainly be interesting to consider the case of more degrees of freedom, from both the dynamical and the control point of view. In principle, all of the dynamical integrity arguments can be repeated word by word, but the actual computation of the safe basin is far from being trivial, because one works in hyperspaces, and it must be faced appropriately. However, the matter of control can still be pursued, as already done -based on a refined theoretical treatment of the homoclinic bifurcation -in the limit case of a continuous (i.e., infinite d.o.f.) system, for the specific case of a buckled beam [START_REF] Lenci | Control of the homoclinic bifurcation in buckled beams: infinite-dimensional vs reduced order modeling[END_REF]. The results of this work call for an extension to more general situations.

Figure 1 .

 1 Figure 1. Helmholtz oscillator: Erosion profiles for the three excitations considered. (a) 4 6 0170, (b) 4 6 0181, (c) 4 6 0185 and (d) 4 6 0190. The vertical segments are the critical thresholds for homoclinic bifurcation.

Figure 2 .

 2 Figure 2. Duffing oscillator: Erosion profiles for 4 6 0180 and (a) harmonic, (b) global optimal with one superharmonic, and (c) right one-side optimal with one superharmonic, excitations. The vertical segments are the critical thresholds for homoclinic bifurcation.

Figure 3 .

 3 Figure 3. Rigid block: erosion profiles of various integrity measures with harmonic excitation and for 4 6 315, 6 0102 and r 6 0195. The vertical line is the critical threshold for homoclinic bifurcation.

Figure 4 .

 4 Figure 4. Rigid block: (a) classical safe basin and (b) "true" safe basin for 2 1 6 0120, 4 6 315, 6 0102 and r 6 0195. The circles for evaluating the Integrity Factor are also reported.

Figure 5 .

 5 Figure 5. Rigid block: erosion curves for 6 015, r 6 017, 4 6 115. The vertical dotted lines are the critical thresholds for homoclinic bifurcation.

Figure 7 .

 7 Figure 7. MEMS: the numerical gains, i.e., the I.F. and G.I.M. normalized to 1 in correspondence with the harmonic excitation, as a function of the superharmonic relative amplitude 2 2 /2 1 . 4 6 017 and 2 1 6 010025.


Jonathan Viquerat

Claire Scheid

Curvilinear DGTD method for nanophotonics applications

Book of abstracts of the th International Conference on Advanced Computational Methods in Engineering, ACOMEN -June .

Curvilinear DGTD method for nanophotonics applications Jonathan Viquerat 1 and Claire Scheid 1,2

1 Nachos project-team, INRIA Sophia Antipolis Méditerranée 2 J.-A. Dieudonné laboratoy, University of Nice Sophia Antipolis e-mails: jonathan.viquerat@inria.fr, claire.scheid@unice.fr

Abstract

Classical nite element methods rely on tessellations composed of straight-edged elements mapped linearly from a reference element, on domains which physical boundaries are indi erently straight or curved. This approximation represents serious hindrance for high-order methods, since they limit the precision of the spatial discretization to second order. Thus, exploiting an enhanced representation of the physical geometry of a considered problem is in agreement with the natural procedure of high-order methods, such as the discontinuous Galerkin method. In the latter framework, we propose and validate an implementation of a high-order mapping for tetrahedra, and then focus on speci c nanophotonics setups to assess the gains of the method in terms of memory and performances.

Key words: Discontinuous Galerkin, curvilinear elements, Maxwell equations, nanophotonics

Curvilinear DGTD formulation

Classical discontinuous Galerkin time-domain (DGTD) methods rely on a linear mapping from a straight-edged reference element to each physical element of the mesh to evaluate the expressions of the nite-element matrices : this allows to save a lot in terms of computational e ciency and memory consumption. Indeed, in the linear case, the nite element matrices for the physical elements are simply multiples of the precalculated matrices of the reference element, since the Jacobian of the corresponding transformation is a constant. In a curvilinear setting, the reference element is mapped to the physical element via a quadratic form, thus allowing a quadratic representation of boundaries. Therefore, the Jacobian of this transformation is no longer a constant, and the matrices have to be evaluated by means of numerical integration, and stored for each physical curved tetrahedron. E cient quadrature and cubature rules can be easily found up to su cient order to our purposes.

A DGTD scheme accouting for curved elements was formulated and implemented in the framework of Maxwell's equations, using centered numerical uxes. A validation step was conducted to verify the stability and accuracy of the method. Realistic situations related to the nanophotonics eld will then be considered that demonstrate the potential of the approach, such as realistically-rounded metallic nanocubes described by a multipole dispersive model.

Acknowledgements

We would like to thank the French DGA (Defence procurement agency) for partial funding of J. Viquerat's PhD thesis.

Validation

A spherical cavity of unit radius is considered, with PEC boundary conditions. A (0, 1, 1) mode is propagated inside it, which exact solution is known. Four di erent rectilinear meshes of increasing re nement were generated in order to check for h-convergence. P 1 to P 4 polynomial approximations were used. For each simulation, the L 2 error is calculated over the whole mesh, and the maximum error levels is retained. The obtained h-convergence results are displayed on table . As can be seen, the use of curvilinear tetrahedra restores quasi-optimal rates. Moreover, it allows to save a lot in terms of degrees of freedom, and therefore in CPU time. Indeed, the curvilinear M solution is three times faster and occupies three times less memory than the linear M solution, for a roughly similar error level. As an example, the P 2 numerical solution on the M and M meshes are displayed on gure .


Stéphane Lanteri
email: stephane.lanteri@inria.fr

Raphaël Léger

Claire Scheid

Jonathan Viquerat

Development of DGTD solver for nanophotonics applications

During the last ten years, the discontinuous Galerkin timedomain (DGTD) method has progressively emerged as a viable alternative to well established finite-difference timedomain (FDTD) and finite-element time-domain (FETD) methods for the numerical simulation of electromagnetic wave propagation problems in the time-domain. We discuss here about the development and application of such a DTGD method for solving the system of time-domain Maxwell equations coupled to material models relevant to nanophotonics. Our efforts aim at improving the accuracy, flexibility and efficiency of the method in view of the numerical treatement of realistic problems.

Generalities about the DGTD method

The DGTD method can be considered as a finite element method where the continuity constraint at an element interface is released. While it keeps almost all the advantages of the finite element method (large spectrum of applications, complex geometries, etc.), the DGTD method has other nice properties which explain the renewed interest it gains in various domains in scientific computing:

-It is naturally adapted to a high order approximation of the unknown field. Moreover, one may increase the degree of the approximation in the whole mesh as easily as for spectral methods but, with a DGTD method, this can also be done very locally. In most cases, the approximation relies on a polynomial interpolation method but the method also offers the flexibility of applying local approximation strategies that best fit to the intrinsic features of the modeled physical phenomena.

-When the discretization in space is coupled to an explicit time integration method, the DG method leads to a block diagonal mass matrix independently of the form of the local approximation (e.g the type of polynomial interpolation). This is a striking difference with classical, continuous finite element formulations. Moreover, the mass matrix is diagonal if an orthogonal basis is chosen.

-It easily handles complex meshes. The grid may be a classical conforming finite element mesh, a nonconforming one or even a hybrid mesh made of various elements (tetrahedra, prisms, hexahedra, etc.).

The DGTD method has been proven to work well with highly locally refined meshes. This property makes the DGTD method more suitable to the design of a hp-adaptive solution strategy (i.e. where the characteristic mesh size h and the interpolation degree p changes locally wherever it is needed).

-It is flexible with regards to the choice of the time stepping scheme. One may combine the discontinuous Galerkin spatial discretization with any global or local explicit time integration scheme, or even implicit, provided the resulting scheme is stable.

-It is naturally adapted to parallel computing. As long as an explicit time integration scheme is used, the DGTD method is easily parallelized. Moreover, the compact nature of methd is in favor of high computation to communication ratio especially when the interpolation order is increased.

As in a classical finite element framework, a discontinuous Galerkin formulation relies on a weak form of the continuous problem at hand. However, due to the discontinuity of the global approximation, this variational formulation has to be defined at the element level. Then, a degree of freedom in the design of a discontinuous Galerkin scheme stems from the approximation of the boundary integral term resulting from the application of an integration by parts to the element-wise variational form. In the spirit of finite volume methods, the approximation of this boundary integral term calls for a numerical flux function which can be based on either a centered scheme or an upwind scheme, or a blend of these two schemes. In the early 2000's, DGTD methods for time-domain electromagnetics have been first proposed by mainly three groups of researchers. One of the most significant contributions is due to Hesthaven and Warburton [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations[END_REF] in the form of a high order nodal DGTD method formulated on unstructured simplicial meshes. The proposed formulation is based on an upwind numerical flux, nodal basis expansions on a triangle (2D case) and a tetrahedron (3D case) and a Runge-Kutta time stepping scheme. In [START_REF] Kabakian | Unstructured gridbased discontinuous Galerkin method for broadband electromagnetic simulations[END_REF], Kakbian et al. describe a rather similar approach. More precisely, the authors develop a parallel, unstructured, high order DGTD method based on simple monomial polynomials for spatial discretization, an upwind numerical flux and a fourth-order Runge-Kutta scheme for time marching. The method has been implemented with hexahedral and tetrahedral meshes.

Finally, a high order nodal DGTD method formulated on unstructured simplicial meshes has also been proposed in the same time frame by Fezoui et al. [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF]. However, contrary to the DGTD methods discussed in [START_REF] Hesthaven | Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell's equations[END_REF] and [START_REF] Kabakian | Unstructured gridbased discontinuous Galerkin method for broadband electromagnetic simulations[END_REF], the method proposed in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] is non-dissipative thanks to a combination of a centered numerical flux with a second-order leap-frog time stepping scheme. The extension of this DGTD method to higher-order leap-frog schemes is studied in [START_REF] Fahs | A high-order non-conforming discontinuous Galerkin method for time-domain electromagnetics[END_REF].

1.1. DGTD method for time-domain nanophotonics

Overview of related works

Numerical modeling of electromagnetic wave propagation in interaction with metallic nanostructures at optical frequencies requires to solve the system of Maxwell equations coupled to appropriate models of physical dispersion in the metal. The most used are the Drude and Drude-Lorentz models. Nevertheless, although Drude and Drude Lorentz dispersion models are widely used among the numerical physicist community and proved their efficiency, they can fail to describe some materials (e.g. transition metals [START_REF] Ulbricht | Carrier dynamics in semiconductors studied with time-resolved terahertz spectroscopy[END_REF]- [START_REF] Wolff | Simple magneto-optic transition metal models for timedomain simulations[END_REF] and graphene [START_REF] De Abajo | Graphene nanophotonics[END_REF]). Furthermore at some scales, nonlocal effects starts to play an important role [START_REF] Moreau | Impact of nonlocal response on metallodielectric multilayers and optical patch antennas[END_REF]. When dealing numerically with Drude and Drude-Lorentz models, the FDTD method is a widely used approach for solving the resulting system of PDEs [START_REF] Taflove | Computational electrodynamics: the finite-difference time-domain method -3rd ed[END_REF]. However, for nanophotonic applications, the space and time scales, in addition to the geometrical characteristics and the physical parameters of the considered nanostructures (or structured layouts of the latter), are particularly challenging for an accurate and efficient application of the FDTD method. Thus, with all their features (as described above), DG methods seem to be well suited in the context of nanophotonics. They are still not prominent compared to FDTD or FETD methods. Nevertheless, recently, unstructured mesh based methods have been developed and have demonstrated their potentialities for being considered as viable alternatives to the FDTD method [START_REF] Ji | High-order DGTD method for dispersive Maxwell's equations and modelling of silver nanowire coupling[END_REF]- [START_REF] Niegemann | Higher-order time-domain methods for the analysis of nano-photonic systems[END_REF]- [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF]- [START_REF] Diehl | Comparison of low-storage Runge-Kutta schemes for discontinuous Galerkin time-domain simulations of Maxwell's equations[END_REF]- [START_REF] Stannigel | Discontinuous Galerkin time-domain computations of metallic nanostructures[END_REF]- [START_REF] Zhu | Analysis of photonic crystals using the hybrid finite-element/finitedifference time domain technique based on the discontinuous Galerkin method[END_REF] for nanophotonics or related applications. Besides, several studies have already been conducted regarding the development of DGTD methods for dispersive media, such as, [START_REF] Lu | Discontinuous Galerkin methods for dispersive media and lossy Maxwell's equations and PML boundary conditions[END_REF]- [START_REF] Lin | Discontinuous Galerkin time-domain method for GPR simulation in dispersive media[END_REF]- [START_REF] Ji | High-order DGTD method for dispersive Maxwell's equations and modelling of silver nanowire coupling[END_REF]- [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF]- [START_REF] Gedney | A discontinuous Galerkin finite element time-domain method modeling of dispersive media[END_REF] to cite a few. Noteworthy, all these studies adopt a DGTD method with upwind numerical fluxes. Furthermore one can find more studies focused on numerical analysis aspects concerning dispersive media [START_REF] Li | Mathematical and numerical study of wave propagation in negative-index material[END_REF]- [START_REF] Li | Numerical convergence and physical fidelity analysis for Maxwell's equations in metamaterials[END_REF].

A non-dissipative DGTD method

Towards the general aim of being able to consider concrete physical situations relevant to nanophotonics, one has to take into account in the numerical treatment, a better description of the propagation of waves in realistic media. The physical phenomenon that on has consider in the first instance here is dispersion. In the presence of an electric field the medium cannot react instantaneously and thus presents an electric polarization of the molecules or electrons that itself influences the electric displacement. In the case of a linear homogeneous isotropic media, there is a linear relation between the applied electric field and the polarization. However, above some range of frequencies (depending on the considered material), the dispersion phenomenon cannot be neglected and the relation between the polarization and the applied electric field becomes complex. In practice, this is modeled by a frequency dependent complex permittivity. Several such models for the characterization of the permittivity exist; they are established by considering the equation of motion of the electrons in the medium and making some simplifications.

There are mainly two ways of handling the frequency dependent permittivity in the framework of time-domain simulations, both starting from models defined in the frequency domain. A first approach is to introduce the polarization vector as an unknown field through an auxiliary differential equation which is derived from the original model in the frequency domain by means of an inverse Fourier transform. This is called the Direct Method or Auxiliary Differential Equation (ADE) formulation. Let us note that while the new equations can be easily added to any time-domain Maxwell solver, the resulting set of differential equations is tied to the particular choice of dispersive model and will never act as a black box able to deal with other models. In the second approach, the electric field displacement is computed from the electric field through a time convolution integral and a given expression of the permittivity which formulation can be changed independently of the rest of the solver. This is called the Recursive Convolution Method (RCM).

We have recently adapted the DGTD-P p method initially introduced in [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF] to deal with various dispersion models. An ADE formulation has been adopted. The resulting ADE-based DGTD method is detailed in [START_REF] Viquerat | Theoretical and numerical analysis of local dispersion models coupled to a discontinuous galerkin time-domain method for Maxwell's equations[END_REF] where we also study the stability and a priori convergence of the method. We first considered the case of Drude and Drude-Lorentz models and, further extend the proposed ADEbased DGTD method to be able to deal with a generalized dispersion model in which we make use of a Padé approximant to fit an experimental permittivity function. The numerical treatment of such a generalized dispersion model is also presented in [START_REF] Viquerat | Theoretical and numerical analysis of local dispersion models coupled to a discontinuous galerkin time-domain method for Maxwell's equations[END_REF]. In this talk, we will discuss about our efforts regarding this DGTD method in order to improve its accuracy, flexibility and efficiency in view of the numerical treatement of large-scale nanophotonics applications. Numerical results will be presented for several of 3D problems ranging from academic test problems to more realistic configurations.


S Lanteri

Claire Scheid

J Viquerat

Discontinuous Galerkin Time-Domain method for nanophotonics

published or not. The documents may come

Drude model

The Drude model describes the response of certain dispersive media to an electromagnetic wave propagating in a certain range of frequencies. The considered model permits to establish a dependency between the permittivity of the material and the angular frequency of the electromagnetic wave in the following form:

ε r,d (ω) = ε ∞ - ω 2 d ω 2 + iωγ d ,
where ω d , γ d and ω are respectively the plasma frequency and the damping constant of the medium, and the angular frequency. Adding a Drude dispersion model therefore implies a coupling, in the time domain, between the electric field E and an additional field, the dipolar current (J p ), through an ODE whose solution we choose is here approximated in a DG framework. A centered fluxes DG method has been chosen to develop a numerical approximation of the problem, given the geometry and the inhomogeneous media to be considered. It is associated with a second-order Leap-Frog scheme in time, therefore inducing a non-dissipative scheme. A theoretical study of the latter has been made, demonstrating an error convergence in O h min(s,p) + ∆t 2 , where p is the spatial order of approximation, and s is related to the regularity assumptions made on the electromagnetic field.

Generalized dispersive model

Recently, several arbitrary dispersive models have been proposed, such as the Critical Points (CP) [START_REF] Vial | A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the FDTD method[END_REF] and the Complex-Conjugate Pole-Residue Pairs (CCPRP) [START_REF] Han | Model dispersive media in FDTD method with complexconjugate pole-residue pairs[END_REF]. Here, another formulation is considered: in accordance with the fundamental theorem of algebra, the permittivity function is written as a decomposition of a constant, one zero-order pole (ZOP), a set of first-order generalized poles (FOGP), and a set of second-order generalized poles (SOGP). This leads to the following expression in the frequency domain:

ε r,g (ω) = ε ∞ - σ jω - l∈L 1 a l jω -b l - l∈L 2 c l -jωd l ω 2 -e l + jωf l .
This general writing allows an important flexibility for two reasons: (i) it unifies most of the common dispersion models in a single formulation (such as Drude, Drude-Lorentz and Debye media for example), and (ii) it permits to fit any experimental data set in a reasonable number of poles (and thus a reasonable number of coefficients). The stability and convergence properties of the resulting DG formulation seem at first glance to be a feasible extension of the Drude case, and are currently under investigation.

Numerical validations

First, a simple validation case has been set up in order to verify that the orders of convergence theoretically obtained for the Drude model based DGTD method were achieved. A unitary PEC cavity with dispersive properties defined by a Drude model is defined. A forced source current has been added to the set of equations in order to obtain an analytical solution to compare with. The calculated convergence orders match the theoretical prediction, as presented in table 1.

Then, a more physical case is considered: a gold nanosphere of radius 20 nm, whose properties are described by a Drude model, is illuminated by a plane wave, the latter being modulated in time by a gaussian function. The discrete Fourier transform of the field is processed along the computation, and is in the end compared with the Mie solution1 of the problem, which is taken as a reference solution. A good adequation has been found between the reference and the numerical solution, which is displayed on figure 1. Further accuracy and performance assessment will be performed on this fondamental case, then paving the way to situations of high interest in the nanophotonic domain.

Figure 1 :

 1 Figure 1: E field magnitude in the frequency domain calculated by the DGTD method in the vicinity of the gold nanosphere. The view is a z = 0 slice.

Table 1 :

 1 Convergence orders with P 1 and P 2 approximations for the PEC cavity case.

The Mie theory provides an analytical solution to diffraction problems in the case of spherical obstacles.