Keywords: Kriging, system inversion, experimental design, sequential design, Tsallis entropy

A design method presented in a previous paper for the sequential generation of observation sites used for the inversion of a prediction model is extended to cope with practical issues such as delayed observations and design of batches of imposed size. The final objective of the construction is to be able to associate with any target T in the output space a value x T of the input factors such that the response of the system at x T will be "close" to T (from an industrial point of view, x T corresponds to manufactory conditions that yield a product whose feature of interest is described by T ). The problem is thus much different from the more standard one where one wishes to build a precise model over the whole input space: here the model only has to be precise over a set of values x T that permit to reach any target T , that is, the observation sites should not be spread over the entire admissible input space, but should rather concentrate in areas that cover the reachable output space when mapped by the system. Examples in low dimensions are presented that illustrate the behavior of the method and allow a comparison to be made with a standard sequential method for designing exploratory experiments.

A Sequential Design Method for the Inversion of an Unknown System

Régis Bettinger * , * * Pascal Duchêne * Luc Pronzato * *

INTRODUCTION

We shall design an experiment (that is, choose observation sites X 1 , . . . , X N ) for a system S : X ⊂ R m → R p , for which no prior information is available. The input domain X is a compact subset of R m and the reachable output domain Y = S(X ) ⊂ R p is unknown. The objective of the design is to construct a model that will allow us to "invert" the system in the following sense: with any given reachable target vector T ∈ S(X ) one wishes to associate an input vector x T such that the value S(x T ) is as close to T as possible. The number N of observations to be made is fixed in advance and the experimental design problem consists in selecting N points X 1 , . . . , X N ∈ X such that the inverse problem above is solved with maximum accuracy after the observation of Y i = Y (X i ), i = 1, . . . , N .

The input/output relationship is modeled by kriging (see, e.g., [START_REF] Sacks | Design and analysis of computer experiments[END_REF]; [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]), which gives, at low computational cost, a prediction and an estimate of the prediction error at each point of the input domain (see Section 2). We use the Matlab toolbox DACE [START_REF] Lophaven | DACE, a Matlab kriging toolbox[END_REF]. After having observed the output values Y N = (Y 1 , . . . , Y N ) ⊤ at the observation sites X 1 , . . . , X N the posterior distribution D(Y (x)|Y N ) of the output at any point x ∈ X is easily constructed. One can then define

x T = arg min x∈X IE{ Y (x) -T 2 |Y N } (1)
= arg min x∈X y N (x) -T 2 + σ 2 N (x) , the inverse prediction for T , with y N (x) = IE{Y (x) |Y N } and σ 2 N (x) = IE{ Y (x)y N (x) 2 |Y N } respectively the model prediction and variance at x. A natural measure of the prediction accuracy for the target T is thus ζ N (T ) = y N (x T ) -T 2 + σ 2 N (x T ) , and two global criteria for the prediction accuracy over the whole reachable output space Y can be considered, J M,N (X 1 , . . . , X N ) = max

T ∈Y ζ N (T ) , J I,N (X 1 , . . . , X N ) = Y ζ N (T ) µ(dT ) ,
with µ(•) some given probability measure on Y indicating the relative interest among target values. A straightforward formulation of the design problem would consist in minimizing IE{J M,N (x 1 , . . . , x N )} or IE{J I,N (x 1 , . . . , x N )} with respect to x 1 , . . . , x N . This is, however, a formidable task and we need to follow another route (the quantities J M,N (X 1 , . . . , X N ) and J I,N (X 1 , . . . , X N ) can nevertheless be used to evaluate a posteriori the quality of a given design in terms of inverse prediction performance, see [START_REF] Bettinger | Design of experiments for response diversity[END_REF]).

Compared to [START_REF] Bettinger | Design of experiments for response diversity[END_REF], we consider a more natural extension to the case when responses are available with delay and several design points must be chosen simultaneously. The idea used in this previous paper consists in choosing observation sites that sample Y as uniformly as possible, in the sense that the N observed responses Y i = Y (X i ), i = 1, . . . , N , are as spread as possible in R p . Since this implies the exploration of the whole domain Y, we expect any possible target T ∈ Y to be "close" to a certain Y i . At the same time, one wishes that the X i 's remain as concentrated as possible in X : indeed, model predictions must be precise for values x T that permit to reach any target T and those inverse predictions x T , obtained by (1), will preferably be chosen in areas where σ 2 N (x) is small, that is, in the neighborhood of some X i 's. This concentration of design points in the input space can be obtained by constructing the design sequentially when the precision of the prediction is properly taken into account at each step, thereby favoring the choice of X n+1 at step n close to points X i already sampled. Choosing the point X n+1 after Y 1 , . . . , Y n have been observed has also the advantage of using all the (increasing) information available and summarized in the current distributions D(Y (x)|Y n ), x ∈ X . An initial design with N 0 points is chosen, with N 0 small compared to N . Since no prior information on the system is available, a space-filling design is used, ensuring that the first N 0 sites are as well spread in X . As suggested by [START_REF] Morris | Exploratory designs for computational experiments[END_REF], we use a latine hypercube design, easy to generate, combined with the optimization of a maximin-distance criterion, see Section 4.2.

Two methods are proposed in [START_REF] Bettinger | Design of experiments for response diversity[END_REF] for choosing X n+1 :

• maximize the conditional (posterior) expectation of the minimum distance between Y (x) and the Y i 's already observed,

X n+1 = arg max x∈X IE{ min i=1,...,n Y (x) -Y i |Y 1 , . . . , Y n } ,
see Section 3.2 (the idea being that, in the average sense, the next output Y n+1 will be "far" from previous ones Y 1 , . . . , Y n ); • maximize the expected second-order Tsallis entropy H 2 [•] [START_REF] Tsallis | Possible generalization of Boltzmann-Gibbs statistics[END_REF] of a kernel density estimator φ n,x formed from the Y i 's and Y (x),

X n+1 = arg max x∈X IE{H 2 [ φ n,x ]} ,
see Section 3.3. The idea here is that maximum entropy over a compact set is obtained for the uniform distribution: we thus expect the method to spread the Y i 's as much as possible in Y (the second-order Tsallis entropy being used because it allows very simple analytical calculations).

Note that the precision of the model prediction is taken into account by both approaches through the conditional distribution of Y (x) used to evaluate the expectation IE{•}.

In [START_REF] Bettinger | Design of experiments for response diversity[END_REF] we showed on a (noise-free) toy example with 2 inputs and 1 output that both approaches yield satisfactory designs.

For practical reasons, it may happen that batches of k points X n+1 , . . . , X n+k must be chosen simultaneously, and that those k points must be selected before the k previous observations are available (knowing, however, the input values X n-k+1 , . . . , X n used in previous batch). This is the case in particular when the analysis of the outputs produced by the system is a time-consuming procedure and external constraints impose that the system S (experimental or simulated) is operated continuously.

A rather crude adaptation of the methods above to such a situation was proposed in [START_REF] Bettinger | Design of experiments for response diversity[END_REF] and illustrated by the application to a 5-inputs/2-outputs noisy system derived from real data collected in oil industry.

We show below that such practical constraints on delayed observations and batch design can easily be taken into account when using the design method based on the maximization of the expected entropy, without requiring any particular approximation (this is not the case, however, for the method based on the maximization of the expected minimum distance). The paper is organized as follows: Section 2 recalls the main properties of kriging predictors. Section 3 presents three sequential design methods, based on the maximization of the kriging variance σ 2 n (x), of the expected minimum distance and of the expected entropy, with their extensions to the "batch-delayed" case. In Section 4 the methods are compared on low-dimensional noise-free toy examples. The adaptation to the case of noisy observations, together with the need of specific optimization algorithms, are mentioned in a conclusion section.

PREDICTION BY KRIGING

For the sake of simplicity of the presentation we only consider the case of a single (scalar) response Y (x). When several responses are present, we thus assume that they are modelled independently. Notice, however, that a more sophisticated approach (co-kriging) permits to take possible correlation between responses into account, see, e.g., [START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty[END_REF]. The model used for Y (x) is

Y (x) = f ⊤ (x)β + Z(x) + ε(x) , where f (x) = (f 1 (x), . . . , f k (x)) ⊤ is a vector of known regression functions, β = (β 1 , . . . , β n ) ⊤ is a vector of unknown parameters, Z(•)
is a second-order stationary stochastic process with zero mean, and ε(x) denotes the observation error at x. Those errors are assumed to be centered, independently distributed with variance σ 2 ε , and independent from Z(x). We note Cov {Z(x 1 ), Z(x 2 )} = σ 2 z ρ(x 1x 2 , ψ) the covariance function of Z, where σ 2 z is the variance and ρ(•) the correlation function of Z, and ψ is a vector of unknown parameters. Popular correlation functions include (with

ψ ∈ R, t ∈ R m ) ρ e (t, ψ) = m i=1 e -ψ|ti| , (2)
ρ g (t, ψ) = e -ψ t 2 , (3)
ρ m (t, ψ) = m i=1 1 Γ(ν)2 ν-1 2 √ ν|t i | ψ ν K ν 2 √ ν|t i | ψ , (4)
called respectively the (isotropic) exponential, Gaussian and Matérn correlation functions, with K ν the modified Bessel function of order ν. The smoothness of the paths of Z is related to the choice of the correlation function: in a nutshell, (2) corresponds to a process with non-smooth trajectories whereas (3) gives infinitely differentiable trajectories. In-between, the regularity of a model built with (4) grows with ν.

Suppose one has collected

n observations Y n = (Y 1 , . . . , Y n ) ⊤ at X 1 , . . . , X n , denote Z n = (Z(X 1 ), . . . , Z(X n ))
⊤ and assume that σ 2 ε and the characteristics of the process ρ(•), σ 2 z and ψ are known. The kriging predictor at an arbitrary point x 0 ∈ X is the best (in the sense of minimum variance) linear unbiased predictor at x 0 . Easy calculations give

y n (x 0 ) = f ⊤ 0 β + r ⊤ 0 R -1 (Y n -F β) ,
(5) where we note

f 0 = f (x 0 ), r 0 = Cor{Z(x 0 ), Z n } ∈ R n , F = (f (X 1 ) . . . f (X n )) ⊤ , and β = (F ⊤ R -1 F ) -1 F ⊤ R -1 Y n , with R = R + (σ 2 ε /σ 2 z ) I n , R = Cor{Z n , Z n } and I n the n-dimensional identity matrix.
The prediction variance (called kriging variance) at x 0 is

σ 2 n (x 0 ) = σ 2 z a ⊤ R a -2a ⊤ r 0 + 1 + σ 2 ε , ( 6
)
where a is obtained from y n (x 0 ) = a ⊤ Y n in (5).

We assume that the process Z(•) is Gaussian and that the observation errors are normal. One can then easily construct the log-likelihood function

l(β, σ 2 z , σ 2 ε , ψ) = - 1 2 n log(σ 2 z ) + log(det R) + (Y n -F β) ⊤ R -1 (Y n -F β) σ 2 z ,
which enable to estimate the unknown kriging parameters β, σ 2 z , σ 2 ε and ψ by maximum likelihood, see, e.g., [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF].

One can easily check that when σ 2 ε = 0 (no measurement errors) y n (X i ) = Y i and σ 2 n (X i ) = 0, i = 1, . . . , n: the predictor is then an interpolator and the prediction variance equals zero at the design points. The method can be put in a Bayesian framework, assuming for instance a prior on β and σ 2 z with σ 2 ε and ψ known. When σ 2 z is also assumed to be known and a non-informative prior is put on β, then the posterior distribution

D[Y 0 |Y n ] of Y 0 = Y (x 0 ) is normal (Y 0 |Y n ) ∼ N ( y n (x 0 ), σ 2 n (x 0 )) .
(7) This is approximately true in other circumstances and can be generalized to the case of simultaneous multiplepoint predictions. Let Y 0 = (Y (x n+1 ), . . . , Y (x n+k )) ⊤ correspond to the vector of responses at the new design sites X 0 = (x n+1 , . . . , x n+k ) and denote Z 0 = (Z(x n+1 ), . . . , Z(x n+k )) ⊤ . Then,

Z n Z 0 ∼ N 0, Σ nn Σ n0 Σ 0n Σ 00
gives (see [START_REF] Pronzato | Robust design with nonparametric models: prediction of second-order characteristics of process variability by kriging[END_REF])

(Y 0 |Y n ) ∼ N ( Y n (X 0 ), Σ n (X 0 )) , (8)
where

Y n (X 0 ) = Σ 0n Σ -1 nn + V ⊤ (F ⊤ Σ -1 nn F ) -1 F ⊤ Σ -1 nn Y n , Σ n (X 0 ) = Σ 00 -Σ 0n Σ -1 nn Σ n0 + V ⊤ (F ⊤ Σ -1 nn F ) -1 V , with V ⊤ = F 0 -Σ 0n Σ -1 nn F , F 0 = (f (x n+1 ), . . . , f (x n+k )) ⊤ .
For k = 1 we recover the normal distribution (7). The distribution (8) will be used to compute expectations when choosing new design points in the presence of batch/delayed observations.

BATCH SEQUENTIAL DESIGN

Maximization of kriging variance

A classical approach for the sequential generation of design points that aim at giving precise kriging predictions over the whole domain X consists in choosing at step n

X n+1 = arg max x∈X σ 2 n (x) , (9)
with σ 2 n (x) given by ( 6). As mentioned above, when there is no observation error we have σ 2 n (X i ) = 0, i = 1, . . . , n, and this approach thus guarantees that no repetition of observations at the same location will take place. (In (9), and also in the optimization problems to follow, the maximization is usually performed by grid search.)

In the case of batch optimization, the choice (9) can be generalized (see [START_REF] Sahama | Sample size considerations and augmentation of computer experiments[END_REF]) into (X n+1 , . . . , X n+k ) = arg max

X0∈X k det Σ n (X 0 ) , (10)
with Σ n (X 0 ) as in ( 8). In the batch-delayed setting, when observations at the k previous sites X n-k+1 , . . . , X n are unknown, we propose to use (X n+1 , . . . , X n+k ) = arg max

X0∈X k det Σ n-k (X 0 ) , (11)
with Σ n-k (X 0 ) the covariance matrix of the random vector (Y (X n-k+1 ), . . . , Y (X n ), Y (x n+1 ), . . . , Y (x n+k )) ⊤ conditional to the observations Y (X 1 ), . . . , Y (X n-k ). Note that when σ 2 ε = 0 this prevents repetitions of observations at the same locations X n-k+1 , . . . , X n as the batch under treatment.

The example used in Section 4 to illustrate the behaviors of the design rules (10, 11) indicates that, although not adapted to the model-inversion problem considered here, those design rules are suitable for generating exploratory space-filling designs in the input space X under the constraint of batch design of imposed size and in presence of delayed observations.

Maximization of minimum distance

Consider first the simplest case k = 1 (where points are added one-at-a-time). The quantity d(y, Y n ) = min i=1,...,n y -Y i measures how far y is from the observations already performed. Using (7), we can compute the expected value of this quantity for an observation made at x and thus choose next design point as

X n+1 = arg max x∈X IE{d(Y (x), Y n ) |Y n } . (12)
An analytic expression for IE{d(Y (x), Y n ) |Y n } is easily obtained when p = 1 ; for larger p, however, the integrant has to be computed numerically, which becomes very timeconsuming.

A generalization to the batch setting (k ≥ 1) could be based on (X n+1 , . . . , X n+k ) = arg max

X0∈X k IE min i=1,...,n+k, j=n+1,...,n+k, j =i Y i -Y j Y n (13)
with Y i = Y (x i ) for i = n + 1, . . . , n + k. However, this is very cumbersome to implement and we do not go any further here; one may refer to [START_REF] Bettinger | Design of experiments for response diversity[END_REF] for a simple heuristic extension of (12) in the batch-delayed setting.

Maximization of Tsallis entropy

This approach follows a similar idea to the method above, but measures the dispersion of the observations in the output space through entropy: since maximum entropy over a compact set is obtained for the uniform distribution, we maximize the (expected) entropy of the distribution of the observations {Y 1 , . . . , Y n , Y (x)} in Y.

First, a kernel-density estimator of this distribution is computed (see, e.g., [START_REF] Wand | Kernel Smoothing[END_REF])

φ n,x (y) = 1 n + 1 n i=1 ϕ Yi,h (y) + ϕ Y (x),h (y) ( 14
)
with ϕ z,h (•) the probability density function (p.d.f.) of the normal N (z, h 2 ). Usually, the smoothing parameter h (window width) is taken decreasing with n in order to ensure good asymptotic properties for the density estimator. In the context considered here, however, the total number of observations N is small, and we keep h constant (h = 0.01, with the range of observed outputs normalized to [0, 1] at each step).

Second, the entropy of ( 14) is computed. A natural candidate is Shannon entropy, which, for φ(•) a p.d.f. on R p , can be written as

H 1 (φ) = -R p φ(t) log[φ(t)] dt.
However, this entropy functional does not yield an analytic formula for a mixture of normal distributions such as ( 14). We thus consider the second-order Tsallis entropy of the p.d.f. φ(•), H 2 (φ) = 1 -R p φ 2 (t) dt, which can be given an analytic expression when substituting the estimator ( 14) for φ. Indeed, we have

H 2 [ φ n,x ] = 1 - 1 (n + 1) 2 n+1 i,j=1 ϕ Yi,h √ 2 (Y j ) ,
where we denote Y n+1 = Y (x). We thus obtain the selection rule

X n+1 = arg max x∈X IE{H 2 [ φ n,x ] |Y n } , (15)
with Y (x) having the distribution (7). One can notice that, moreover, IE{H 2 [ φ n,x ] |Y n } can be given an analytic expression using the property R p ϕ a,σ (t) ϕ b,δ (t) dt = ϕ a, √ σ 2 +δ 2 (b). In the batch setting, we replace the kernel estimator ( 14) by

φ n,X0 (y) = 1 n + k n i=1 ϕ Yi,h (y) + n+k i=n+1 ϕ Y (xi),h (y)
and obtain the multiple-point selection rule (X n+1 , . . . , X n+k ) = arg max

X0∈X k IE H 2 φ n,X0 Y n ,
(16) which can still be given an analytic expression. In the batch-delayed case, we use an approach similar to that in Section 3.1 and choose (X n+1 , . . . , X n+k ) as (X n+1 , . . . , X n+k ) = arg max

X0∈X k IE H 2 φ n,X0 Y n-k . (17)

EXAMPLES

We first consider a 1-input/1-output example to illustrate graphically the behaviors of selection rules (16, 17) and (10, 11) and compare their performances. We then consider the same 2-inputs/1-output example as in [START_REF] Bettinger | Design of experiments for response diversity[END_REF] and generate the design with the rule ( 16), which allows some comparison to be made with the case when the points are added one-at-a-time with (15).

The kriging parameters are estimated by maximum likelihood, excepted for ν which is set to 3.7.

A toy example with

m = p = 1
We illustrate the behaviors of the selection rules on a given set of noise-free scalar observations depending on a scalar input variable x ∈ X , with X formed of 101 equally spaced points in [0, 1].

We first compare the influence of the correlation function on the new set of points proposed by the design rule ( 16) for batches of size k = 2.

The locations of the points proposed by the algorithm is plotted on Figure 1. The n = 6 current observations (set to arbitrary values) are represented by dots, the solid lines correspond to the kriging predictions and the dashed lines indicate the 95% confidence intervals around the predictions. The locations of the two design points generated by ( 16) are indicated by stars. In the exponential case, the model is not smooth and the confidence bounds are very loose; the new points proposed by the design rule are close to previous points, in regions where predictions are reliable. The fact that confidence bounds are unreasonably tight for the Gaussian correlation is known in the kriging litterature, see [START_REF] Stein | Interpolation of Spatial Data[END_REF], and might produce an overconfidence in the predictions, although this is not the case here. Also, Gaussian correlation matrices are often poorly conditioned [START_REF] Ababou | On the condition number of covariance matrices arising in kriging, estimation, and simulation of random fields[END_REF]; since inverses of correlation matrices are required this raises numerical difficulties (even if Cholesky factorization is used). Therefore, in the rest of the paper we use (and in general recommend the use of) the Matérn correlation function (4).

Consider now the behavior of the design rule (10) on the same data set (with the Matérn correlation), see Figure 2. Similarly to its well-known version (9) for k = 1, the design rule (10) enforces exploration and proposes points located in low-sampled areas in the input space. It is therefore not especially adapted to the model-inversion problem considered in this paper.

In Figure 3 we compare the rules ( 17) and ( 11) with the same set of observations as above (dots) and considering the previous optimal batches of size 2 as being under treatment (circles). One may notice that the entropy-based rule (17) (top) places the next two points in areas where predictions are precise and reasonably far from previous observations (or expected observations for the two delayed ones). The variance-based rule (11) (bottom) has the same tendency as the rule (10) to generate new points far from previous design points, thus enforcing a space-filling property in the input domain. Although not adapted to the context considered here, this rule can be of interest in a batch-delayed setting when the objective is to build an accurate model over the entire input domain.

A toy example with m = 2 and p = 1

The system is now supposed to follow the equation (no observation error, σ 2 ε = 0)

Y (x 1 , x 2 ) = 0.2e x1-3 + 2.2|x 2 | + 1.3x 6 2 -2x 2 2 -0.5x 4 2 - 0.5x 4 2 + 2.5x 2 1 + 0.7x 3 1 + sin(5x 1 ) cos(3x 2 1 ) + 3 (8x 1 -2) 2 + (5x 2 -3) 2 + 1 .
The initial design is a 9-points maximin latin hypercube (several such initial designs have been used, showing little influence on the performance), points are then added sequentially by pairs of k = 2 up to N = 19 with the rule ( 16). The input domain corresponds to a regular grid of 31×31 points in the square [0, 1] 2 . (This small number of points permits to obtain the optimal solution at each step 16) in a 2-inputs/1output example: inputs (top) and outputs (bottom).

CONCLUSIONS

This paper proposed extensions of two sequential design rules to situations where runs should be made several-ata-time and practical considerations impose a delay in the availability of observations. Simple examples illustrated that the first one is appropriate when the objective is to construct an accurate model over the whole input domain, while the second one is especially suited for model-inversion.

Both methods guarantee that there will be no repetition of observations at the same design site when there is no measurement error (σ 2 ε = 0). The situation is different when σ 2 ε = 0. In that case, we suggested in [START_REF] Bettinger | Design of experiments for response diversity[END_REF] to use the re-interpolation technique proposed by [START_REF] Forrester | The design and analysis of 'noisy' computer experiments[END_REF] for computer experiments. A first kriging predictor is constructed, taking the presence of measurement errors into account; it does not interpolate the data, that is, y n (X i ) = Y i . A second kriging predictor is then constructed, assuming that there are no observation errors and using the predictions y n (X i ) as if they were observations. One can show that for any x the prediction at this second stage coincides with y n (x) of the first stage, with the noticeable difference that the prediction variance for the second predictor is now zero at the observation points X i , which permits to avoid repetitions.

For both approaches, generating k points at a time requires the solution of a k × m-dimensional optimization problem at each step, which quickly becomes computationally unfeasible when k increases (for the low dimensional examples considered, with m = 1, 2 and k = 2, the input domain corresponded to a regular grid and all possible choices of design points could be tested at each step). The development of suitable algorithms (e.g., of the exchange type) is under current investigation and forms a prerequisite to the solution of the 5-inputs/2-outputs model-inversion problem presented in [START_REF] Bettinger | Design of experiments for response diversity[END_REF].

Fig. 1 .

 1 Fig. 1. Effect of the correlation function (exponential:top; Gaussian: middle; Matérn: bottom) on the choice of k = 2 new points with the rule (16).

Fig. 4 .

 4 Fig. 4. Design points generated by (16) in a 2-inputs/1output example: inputs (top) and outputs (bottom).

  Figure4illustrates the way points are added at each step. The input values are plotted on the top, with dots representing the initial design, numbers indicating the sequence of batches generated and dashed lines corresponding to contour plots of the response Y (•, •). The associated output values are plotted on the bottom. Note that the points generated are concentrated in the input space; the dispersion of the outputs in Y = [-0.67, 4.75] (bottom) is a bit worse than for k = 1, see[START_REF] Bettinger | Design of experiments for response diversity[END_REF].

	1
	0.8
	0.6
	0.4
	y
	0.2
	0
	-0.2
		design points
	-0.4	predictions 95% conf. int.
		new points
	0	0.2	0.4	x	0.6	0.8	1
	Fig. 2. Choice of k = 2 new points with the rule (10)
	(Matérn correlation).
	1
	0.8
	0.6
	0.4
	y
	0.2
	0
	-0.2	design points
		predictions
	-0.4	95% conf. int. points under way new points	x
	0	0.2	0.4		0.6	0.8	1
	1
	0.8
	0.6
	0.4
	y
	0.2
	0
	-0.2	design points
		predictions
	-0.4	95% conf. int. points under way
		new points
	0	0.2	0.4	x	0.6	0.8	1
	Fig. 3. Choice of k = 2 new points with the rules (17) (top)
	and (11) (bottom) in presence of delayed observations.
	by exhaustive search in a reasonable computing time. We
	noticed, however, that the results are rather sensitive to
	the grid size. Considering finer grids, with an optimization
	algorithm used at each step instead of exhaustive search,
	would thus be of interest.)


David Andersen
email: dandersn@illinois.edu

Yuli Starosvetsky

Alexander Vakakis

Lawrence Bergman

Dynamic instabilities in coupled oscillators induced by geometrically nonlinear damping

Keywords: Nonlinear damping, Essential nonlinearity, Nonlinear energy sink, Impulsive orbit

Introduction

Effective suppression of unwanted vibrations is a topic of special interest in various engineering applications. One of the well-known solutions for this problem is a light-weight linear attachment, coupled to an externally excited linear or weakly nonlinear substructure via viscous damping and a linear spring [START_REF] Hartog | Mechanical Vibrations[END_REF]. This linear light mass attachment is generally referred in the literature as a tuned mass damper. The purpose of this attachment is to attenuate vibrations over a narrow frequency range centered at the natural frequency of the absorber. The effective bandwidth is governed by the damping in the absorber, and a trade-off exists between attenuation efficiency and bandwidth. De-spite the aforementioned advantages of the linear absorber, it has a major drawback: namely, it is only effective in the neighborhood of a single frequency. Recently, it was demonstrated that the use of strongly nonlinear attachments can provide a tool for achieving dynamical regimes that are unavailable in common weakly nonlinear systems and also provide fairly efficient broadband suppression unlike the tuned mass damper. The conceptually new phenomenon of targeted energy transfer was observed in these strongly nonlinear systems, where either a linear or weakly nonlinear substructure is coupled via a strongly nonlinear (nonlinearizable) spring and viscous damping with a lightweight mass attachment. The light mass attachment and the nonlinear coupling are collectively referred to in the literature as a nonlinear energy sink (NES). The observed phenomenon of targeted energy transfer, where a substantial amount of energy initially stored in the linear substructure is transferred to the NES in a nearly one-way irreversible fashion, was intensively studied in recent publications [START_REF] Gendelman | Transition of energy to nonlinear localized mode in highly asymmetric system of nonlinear oscillators[END_REF][START_REF] Gendelman | Energy pumping in nonlinear mechanical oscillators. I. Dynamics of the underlying Hamiltonian system[END_REF][START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators[END_REF][START_REF] Vakakis | Inducing passive nonlinear energy sinks in linear vibrating systems[END_REF][START_REF] Vakakis | Dynamics of linear discrete systems connected to local essentially nonlinear attachments[END_REF][START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF]. Although current observations of targeted energy transfer have been limited to certain types and magnitudes of excitation, this phenomenon shows great promise for future engineering applications. It is interesting to note that recently there has been a growing interest toward the applications of NESs for the purposes of passive sound control [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF] and vibration isolation [START_REF] Gatti | On the response of a harmonically excited two degree-of-freedom system consisting of a linear and a nonlinear quasi-zero stiffness oscillator[END_REF][START_REF] Gatti | On the interaction of the responses at the resonance frequencies of a nonlinear two degrees-of-freedom system[END_REF].

There have been several recent studies of the dynamics of an impulsively-excited linear substructure with a NES attached; however, these restricted the nonlinear coupling term solely to the stiffness element [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF][START_REF] Lee | Triggering mechanisms of limit cycle oscillations due to aeroelastic instability[END_REF][START_REF] Kerschen | Irreversible passive energy transfer in coupled oscillators with essential nonlinearity[END_REF]. Even with only nonlinear stiffness coupling, the system dynamics were very complicated. These works identified that particular impulse loadings resulted in a large amount of unidirectional targeted energy transfer from the linear oscillator (LO) to the NES. Also introduced by Lee et al. [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF][START_REF] Lee | Triggering mechanisms of limit cycle oscillations due to aeroelastic instability[END_REF] and Kerschen et al. [START_REF] Kerschen | Irreversible passive energy transfer in coupled oscillators with essential nonlinearity[END_REF] was the concept of the frequencyenergy plot (FEP), which was used to determine the types of resonance captures that the dynamics can support, and also provide global insight into the transient resonant transitions between various nonlinear normal modes of the system for decreasing energy due to damping.

It is important to note that all the previous studies of the dynamics of an impulsively-excited, two degreeof-freedom system considered linear viscous damping in the coupling between the linear oscillator and the NES. In this paper, we will be examining the dynamics of the system with both nonlinear damping and stiffness coupling elements in order to show that the presence of nonlinear damping can introduce complex dynamical phenomena, including transient instabilities. Moreover, these strongly nonlinear damping terms are physically realizable through the geometry of the deformation of linear elements. The physical realization of nonlinear stiffness was examined in Gourdon et al. [START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF] and Gatti et al. [START_REF] Gatti | On the response of a harmonically excited two degree-of-freedom system consisting of a linear and a nonlinear quasi-zero stiffness oscillator[END_REF][START_REF] Gatti | On the interaction of the responses at the resonance frequencies of a nonlinear two degrees-of-freedom system[END_REF] and the concept of nonlinear damping and its physical realization was first introduced by Koumousis and Kefala [START_REF] Koumousis | On the dynamic behavior of a lightweight isolator for museum artifacts[END_REF]. A preliminary study of the dynamic effects of nonlinear damping is given in the works by Triplett et al. [START_REF] Triplett | Energy harvesting from an impulsive load with essential nonlinearities[END_REF] and Quinn et al. [START_REF] Quinn | Comparing linear and essentially nonlinear vibration-based energy harvesting[END_REF].

The structure of this paper is as follows: Sect. 2 is devoted to the model formulation and also describes the physical realization of geometrically nonlinear damping in the system under investigation. The dynamics of the underlying Hamiltonian system (with zero damping) and the construction of the FEP are discussed in Sect. 3, whereas in Sect. [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators[END_REF] we study damped transitions of the system with nonlinear damping by superimposing the corresponding wavelet spectra to the Hamiltonian FEP. The numerical results provide ample motivation for the theoretical study of the fundamental resonances corresponding to 1:1 and 1:3 resonant conditions carried out in Sects. 5 and 6,r espectively. It is shown that, for the case of 1:3 resonance capture, the system under investigation may be approximated by a slow flow equation of lower dimensionality than was used in previous studies; good correspondence of the analytical results with direct numerical simulations is demonstrated. Finally, in Sect. 7, the basic findings of this work are summarized.

Problem formulation

We examine the dynamics of a two degree-of-freedom system of coupled oscillators depicted in Fig. 1.I t is composed of a linear damped oscillator (the lower mass equal to unity) coupled to a light oscillator (the upper mass equal to ǫ) through a viscous damper (λ 2 ) and an additional configuration of side pairs of springs (C/2) in parallel to viscous dampers (λ 3 /2). In this paper, we will work with nondimensional parameters, and thus the results discussed herein are realizable through the proper selection and combination of Fig. 1 Two degree-of-freedom system with essentially nonlinear stiffness and strongly nonlinear damping the physical elements. Even though the side springdamper pairs possess linear force-response characteristics (i.e., the force produced by the side spring is proportional to its axial deformation and the force produced by the side viscous damper is proportional to the rate of its axial deformation), it can be shown [START_REF] Koumousis | On the dynamic behavior of a lightweight isolator for museum artifacts[END_REF][START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF][START_REF] Mcfarland | Experimental investigation of targeted energy transfer in strongly and nonlinearly coupled oscillators[END_REF] that due to the geometry of their configuration the resulting relations between the vertical components of the net forces exerted to the upper mass by these elements are strongly nonlinear functions of the relative vertical response between the two masses. Indeed, the equations governing the vertical oscillations of the system of Fig. 1 are expressed as

x ′′ + λ 1 x ′ + λ 2 (x ′ -v ′ ) + λ 3 (x -v) 2 (x ′ -v ′ ) + ω 2 o x + C(x -v) 3 = F(t) ǫv ′′ + λ 2 (v ′ -x ′ ) + λ 3 (v -x) 2 (v ′ -x ′ ) + C(v -x) 3 = 0 (1)
Hence, the two springs produce an essentially nonlinear (nonlinearizable) stiffness nonlinearity of the third degree, whereas the two damping elements produce an essential nonlinearity of the third degree of the product of the relative velocity times relative response squared. The geometrically nonlinear terms are derived in the Appendix to this report and are approximate; they are valid for large deflections by the selection of a suitably long length L in the physical setting. We conclude that due to the geometry of the motion, the dynamics of the system of Fig. 1 that is constructed out of linear elements, is strongly nonlinear.

Clearly, when the linear damping element in the coupling vanishes, λ 2 = 0, the system possesses purely nonlinear coupling; in the following analysis, we will assume that this system is purely passive by letting λ 3 ≥ 0. Moreover, in the numerical simulations considered in this work the dynamics of system (1) will be analyzed subject to an impulsive excitation F(t)= I o δ(t) on the linear oscillator with the system being initially at rest. This is equivalent to imposing the initial conditions and excitation

x(0) = v(0) = v ′ (0) = 0,x ′ (0) = I o , F(t)= 0 (2)
Finally, we wish to study the dynamics in the limit of a light-weight attachment, so we impose the restriction 0 <ǫ≪ 1 and consider ǫ as the small parameter of the problem; in fact, ǫ will be used as the perturbation parameter in the asymptotic analyses of Sects. 5 and 6.

In order to demonstrate the existence of transient dynamic instabilities induced by the geometrically nonlinear damping, the system described by ( 1) and (2) first will be examined numerically. In Sect. 4, we will depict the results of the numerical simulations by plotting the wavelet transform spectra of the computed time series on the frequency-energy plot (FEP) of the underlying Hamiltonian system (corresponding to zero damping). As discussed in detail in [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF] and in works referenced therein, superimposing wavelet spectra of damped transitions to the plot of periodic (and a subset of quasiperiodic) orbits of the underlying Hamiltonian system can reveal the dominant transient resonance interactions (resonance captures) that occur in the damped dynamics and, therefore, can provide useful insight into the transient behavior of the damped system. Hence, in Sect. 3,w e provide a brief aside on the Hamiltonian FEP of the undamped version of system (1). The numerical simulations of Sect. 4 will provide ample motivation for the analytical work that follows. In the analysis of Sects. 5 and 6, the system of equations will be reduced to the equations governing the slow flow dynamics, and the resulting reduced system will be used to predict and interpret the numerically observed behaviors. Two particular responses will be examined: namely, (a) the bifurcation of the dynamics from a 1:1 to 1:3 resonance capture; and (b) the case where the dynamics tracks the impulsive resonance branch on the FEP followed by a dynamic instability developing as the dynamics is captured into 1:3 resonance.

Review of the dynamics of the Hamiltonian system

First, we consider the dynamics of the Hamiltonian system that results when λ 1 = λ 2 = λ 3 = 0. This system has been described in detail in [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF], and the structure of its periodic and quasiperiodic orbits was systematically studied. This system possesses a very complicated structure of periodic orbits which can be conveniently represented in a frequency-energy plot (FEP) where the frequency index of a periodic orbit is depicted as function of its (conserved) energy. In Fig. 2, we depict the Hamiltonian FEP for parameters ǫ = 0.05, ω o = 1.0, and C = 1.0, which was constructed by solving numerically an appropriately defined nonlinear boundary value problem [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF]. Symmetric periodic orbits Snm± are curves in the configuration plane (x, v), with the ratio (m:n)i ndicating the internal resonance realized, that is, the ratio of the basic frequency of the linear oscillator over that of the nonlinear attachment is equal to m/n. For example, a 1:1 internal resonance is realized on S11±, with both the linear oscillator and the nonlinear attachment oscillating with identical frequencies. The (±) signs indicate in-phase or out-of-phase mo-Fig. 2 Frequency-energy plot (FEP) depicting the periodic orbits of the Hamiltonian system (1); impulsive orbits (IOs) are denoted by bullets (")[(+) are bifurcation points with four Floquet multipliers equal to +1, and (!)with two Floquet multipliers equal to +1andtwoto-1] [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF] tions of the LO and the nonlinear attachment. Unsymmetric periodic orbits Upq± are represented by Lissajous curves in the configuration plane. Moreover, the frequency indices used to characterize periodic orbits on branches Snm± and Unm± are defined by frequency index = nω o /m; the only exceptions are S11± which form the main backbones of the FEP. Certain periodic orbits, termed impulsive orbits-IOs, are of particular importance since they satisfy the special initial conditions x(0) = v(0) = v ′ (0) = 0, and x ′ (0) = X. Hence, an IO represents the response of the Hamiltonian system (initially at rest) after application of an impulse Xδ(t) to the linear oscillator. Three basic types of periodic orbits are discerned in the FEP of Fig. 2.T h ebackbone branches S11± are defined over broad frequency and energy ranges. The subharmonic tongue branches Snm± and Unm± are defined for fixed rational frequency indices and over finite energy ranges. Typically the subharmonic tongues occur in stable-unstable pairs [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF], for example, there are two subharmonic tongues corresponding to 1:3 internal resonance between the LO and the nonlinear attachment, with the in-phase tongue S13+ being unstable and the out-of-phase tongue S13being stable. Finally, the manifold of impulsive orbits is composed of the previously described time-periodic IOs corresponding to rational frequency indices as well as quasiperiodic IOs corresponding to irrational frequency indices. The essential role of the IOs in targeted energy transfer [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF][START_REF] Kerschen | Irreversible passive energy transfer in coupled oscillators with essential nonlinearity[END_REF]insystem(1) will be discussed below. The complexity of the dynamics of the underlying Hamiltonian system (1) is in marked contrast to the dynamics of the limiting linear system obtained as ǫ → 0 (i.e., in the limit of zero mass for the nonlinear attachment), for which the FEP consists of a single horizontal line with frequency index equal to ω o = 1.0; this demonstrates the drastic global effect that the lightweight, local nonlinear attachment has on the Hamiltonian dynamics. Finally, we note that nearly horizontal branches in the FEP correspond to periodic orbits whose frequency depends weakly on energy, so they represent weakly nonlinear motions; on the contrary, curved branches correspond to strongly nonlinear motions since they correspond to periodic orbits whose frequency depends strongly on energy. For example, considering the topology of the branch S11on the FEP, we note that the corresponding out-ofphase periodic orbits satisfying a 1:1 resonance can be either weakly nonlinear (corresponding to the near horizontal sub-branch of S11-) or strongly nonlinear (corresponding to the curved sub-branches of S11-). In the next section, we study the impulsive response of the damped system (1) by superimposing the wavelet transform spectrum of the damped time series to the FEP of Fig. 2 in order to identify the dominant resonance captures occurring in the dynamics.

Damped responses on the Hamiltonian FEP

We now consider the impulsive response of system ( 1), ( 2) when the dissipative terms are included. As discussed in [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF] a way to interpret complex nonlinear transitions in the transient dynamics of dissipative coupled oscillators is to superimpose the wavelet transform spectra of the corresponding time series to the frequency-energy plot (FEP) of the underlying Hamiltonian system. This representation reveals the basic resonance interactions that take place in a given transition [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF].

As shown below the strongly nonlinear dissipative system (1), ( 2) possesses complex dynamics when forced in certain frequency and energy ranges. Unless otherwise noted, in the following simulations we consider the fixed parameters ω o = 2.9, λ 1 = 0.0015, λ 2 = 0, λ 3 = 0.010, C = 1.0, and ǫ = 0.05, so the coupling term between the LO and the attachment possesses purely nonlinear damping. Moreover, we consider the damped responses of system (1), (2) for different strengths of applied impulse, in order to study the effect on the dynamics of the initial energy applied impulsively to the LO.

Starting from Fig. 3 which corresponds to an applied impulse I o = 1.0, we deduce the presence of higher frequency harmonics in the initial, highly energetic phase of the motion; these result due to nonlinear resonance interactions between the attachment and the LO at higher frequency subharmonic tongues [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF] and are similar to what was previously observed in coupled oscillators with analogous configurations but with linear damping. Following the initial high frequency resonance interactions, the damped dynamics follows the branch S11+ (i.e., the dynamics lock into a state of 1:1 resonance) as energy decreases due to the action of the damping elements. This dynamics is similar to that of the system with weak linear damping, so we conclude that in this case no new interesting dynamics results by the addition of the geometrically nonlinear damping. tonian FEP, which indicates a continuous fluctuation of the dynamics between two states of (in-phase) 1:2 and (out-of-phase) 1:3 resonance capture (or equivalently between two states where the LO oscillates two or three times faster than the strongly nonlinear attach-ment). After this period of transient instability, the system settles into a state of 1:3 resonance capture as the dynamics tracks the out-of-phase subharmonic tongue S13-and then decays to zero with increasing time. We note that this instability is attributed solely to the Fig. 5 Impulsive response of the system with linear damping for I o = 0.85: Wavelet spectrum of relative response on the FEP indicating absence of dynamic instabilities passive nonlinear damping of the system, as inferred from the wavelet plot of Fig. 5 corresponding to the system with only linear damping and nonlinear stiffness (i.e., with λ 1 = 0.0015, λ 2 = 0.010, λ 3 = 0, and I o = 0.85 and all other parameters fixed), possessing weakly nonlinear damped dynamics that occur close to the natural frequency of the LO on the linearized branch S11-. In the latter case, the responses of both the LO and the attachment follow the expected (linearized) exponentially damped decays, with complete absence of strongly nonlinear resonance interactions between them. Although dynamic instabilities have been observed for the system of solely nonlinear stiffness, it is clear from Figs. 4 and5 that the nonlinear damper in this case causes the dynamic instability and alters the dynamical behavior of the system.

Further decreasing the strength of the applied impulse we note a strengthening of the transient dynamic instability in the response (cf. Figs. 6 and7). Considering first the response depicted in Fig. 6, we note the presence of the basic harmonic ω o (the natural frequency of the LO) and of a multitude of lower harmonics in the wavelet spectrum of the relative response v(t) -x(t). It is clear that multiple transient dynamic instabilities appear in the response, and these are associated with "bursts" in the frequency content of the transient response at different time periods. We note that during this complex transition the LO oscillates predominantly at its own natural frequency ω o ,sothe rich lower frequency content is mainly due to the oscillation of the strongly nonlinear attachment.

Of particular interest is the initial transient response of the strongly nonlinear attachment (phase I in Fig. 6) with continuously decreasing frequency as its dynamics tracks approximately the manifold of impulsive orbits (IOs) of the underlying Hamiltonian FEP.Atthe same time, the LO oscillates at its own natural frequency. This is followed by a brief out-of-phase 1:3 resonance capture on S13-(phase II), with the LO oscillating at frequency ω o and the strongly nonlinear attachment at ω o /3. Then a sudden transient dynamic instability follows, resulting in a high-frequency "burst" in the attachment response as evidenced in the wavelet spectrum (phase III), and leading to an extended period of in-phase 1:2 resonance capture on branch S12+ (phase IV) during which the LO oscillates twice as fast as the attachment. This is followed by escape from resonance capture and transition of the dynamics to a series of higher-order but briefer resonance captures, with a final lower-frequency dynamic instability (albeit of lesser magnitude) occurring in the process (phase V). We emphasize that, as in the previous case, Fig. 6 Impulsive response of the system with nonlinear damping for I o = 0.70: (a) Relative response v(t) -x(t), (b) wavelet spectrum of relative response, (c) wavelet spectrum on the FEP the complex dynamics is solely due to the presence of nonlinear damping in the system, since no such dynamic phenomena are realized in the corresponding system with linear damping. Moreover, we note the complete absence of 1:1 resonance capture in the dynamics for this particular value of applied impulse.

A stronger dynamic instability induced by the geometric damping nonlinearity is depicted in Fig. 7, corresponding to a further decrease of the magnitude of the impulse to I o = 0.65. In this case, we note three distinct phases of the damped transition. In phase I (cf. Figs. 7a,7b) the dynamics of the strongly nonlinear at- tachment follows closely the manifold of IOs (just as in the previous case), and its frequency decreases continuously with time; during this phase, the LO possesses a single dominant frequency equal to ω o . When the frequency of the attachment reaches the neighborhood of ω o /3, a strong dynamic instability occurs (phase II) as out-of-phase 1:3 resonance capture occurs, and the oscillation of the attachment makes a sudden transition to a higher amplitude. This is followed by escape from resonance capture (phase III) and gradual decay of the amplitude of oscillation of both the LO and the strongly nonlinear attachment. In addition, as in the previous case (but contrary to the cases corresponding to stronger applied impulses examined previously), there is complete absence of 1:1 resonance capture in the dynamics. This particular instability will be examined in detail in Sect. 6.

From the previous results, we conclude that geometrically nonlinear damping gives rise to different types of resonance captures that result in transient dynamic instabilities in the response of system (1), (2); moreover, these instabilities have no counterparts in the corresponding system with weak linear damping. Motivated by these findings, in the following two sec-tions, we provide an analysis of two specific types of resonance captures that can occur in the transient damped dynamics in an attempt to gain an understanding of the mechanisms that govern dynamic instabilities in the system under consideration. By this analysis, we will gain the capacity to predict dynamic instabilities due to geometrically nonlinear damping, or alternatively to use such instabilities in a constructive way in engineering applications, e.g., for shock mitigation or energy harvesting.

Capture into 1:1 resonance

For the analytical study carried out in this and the next sections, we rescale system (1) to the parameter values ω o = 1.0, ǫ = 0.05, λ 1 = λ * 1 ǫ, λ 3 = λ * 3 ǫ, and C = C * ǫ, where the damping and nonlinear stiffness terms are assumed to be of O(ǫ). Moreover, we again consider impulsive excitation of the linear oscillator with the system being initially at rest. We note that the rescalings do not restrict the generality of the analysis. In this section, we wish to study 1:1 resonance captures in the damped dynamics similar to the ones depicted in the previous numerical simulations (see Figs. 3 and4).

For the rescaled parameters and applied impulsive excitation I o = 0.385, the time series of the responses of the LO and the strongly nonlinear attachment for 1:1 resonance capture are shown in Figs. 8a and8b, respectively, and the corresponding wavelet spectrum of the relative response superimposed to the underlying Hamiltonian FEP is presented in Fig. 8c. We note that the damped dynamics mainly follows the branch S11+ (that is, in-phase transient 1:1 resonance capture), before escaping from resonance and making a transition to branch S13 and 1:3 resonance capture. These resonance captures as well as the bifurcation (jump) between them can be examined analytically by reducing the damped dynamics to slow flow equations.

As a matter of convenience, we will introduce the coordinate transformations

u = x + ǫv w = x -v (3)
where u is related to the motion of the center of mass and w to the relative motion between the LO and strongly nonlinear attachment of the system of Fig. 1. Substitution of (3)i n t o( 1) results in the transformed equations of motion

u ′′ + (λ 1 u ′ + u)/(1 + ǫ) + ǫ(λ 1 w ′ + w)/(1 + ǫ) = 0 w ′′ + (λ 1 u ′ + u)/(1 + ǫ) + ǫ(λ 1 w ′ + w)/(1 + ǫ) + λ 3 w 2 w ′ (1 + ǫ)/ǫ + Cw 3 (1 + ǫ)/ǫ = 0 (4)
Next, the complexification-averaging technique, first introduced by Manevitch [START_REF] Manevitch | Complex representation of dynamics of coupled oscillators[END_REF] (and extensively applied to the dynamics of nonlinear coupled oscillators in [START_REF] Vakakis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF]), is applied to system (4) in order to reduce the transient dynamics to an appropriate slow flow. This is the dynamical system governing the ("slow") modulations of the "fast" frequency components and, in essence, is a set of modulation equations. To apply complexification-averaging to study 1:1 resonance capture in this problem, it is necessary to formulate the appropriate ansatz for the transient dynamics; hence, we assume that the LO and the strongly nonlinear attachment are in 1:1 resonance through the entire damped transition possessing identical fast frequencies; moreover, taking into account the form of the transformed equations (4), this common fast frequency must necessarily be equal to unity (the linearized natural frequency of the u-oscillator). To this end, complex functions ψ 1 (t) and ψ 2 (t) are introduced, and are further expressed as

ψ 1 (t) ≡ u ′ + iu = φ 1 (t) exp(it) ψ 2 (t) ≡ w ′ + iw = φ 2 (t) exp(it) (5)
where the φ i (t) are slow complex amplitudes modulating the fast oscillations e it , and i = (-1) 1/2 . Substitution of ( 5)into(4) and subsequent averaging over the fast frequency (equal to unity), we derive the slow flow governing 1:1 resonance capture as

φ ′ 1 + iǫ + λ 1 2(ǫ + 1) φ 1 + ǫ(λ 1 -i) 2(ǫ + 1) φ 2 = 0 φ ′ 2 + 1 + ǫ 8ǫ (λ 3 -3Ci)|φ 2 | 2 φ 2 + λ 1 -i 2(ǫ + 1) φ 1 + i + ǫλ 1 2(ǫ + 1) φ 2 = 0 (6)
The dynamics of the slow flow is now studied by a multiple scales analysis. We introduce the new time scales, τ 0 = t, τ 1 = ǫt,..., and the complex amplitudes in [START_REF] Vakakis | Dynamics of linear discrete systems connected to local essentially nonlinear attachments[END_REF] are expressed in terms of these new time scales which are treated as independent from each other and expressed in power series of ǫ so that φ 1 (τ 0 ,τ 1 ) = φ 10 (τ 0 ,τ 1 ) + ǫφ 11 (τ 0 ,τ 1 )

+••• φ 2 (τ 0 ,τ 1 ) = φ 20 (τ 0 ,τ 1 ) + ǫφ 21 (τ 0 ,τ 1 ) +••• (7)
Expressing the time derivatives in [START_REF] Vakakis | Dynamics of linear discrete systems connected to local essentially nonlinear attachments[END_REF] in terms of derivatives of the new time scales, d dt = ∂ ∂τ 0 + ǫ ∂ ∂τ 1 + •••, substituting (7)i n t o( 6), and matching terms of t h es a m ep o w e r so fǫ, we obtain a hierarchy of subproblems that govern the dynamics of 1:1 resonance capture at different orders of approximation. Considering only O(1) terms in (6), we derive the leading order subproblem

∂φ 10 ∂τ 0 = 0 ⇒ φ 10 = φ 10 (τ 1 ) ∂φ 20 ∂τ 0 + λ * 3 -3C * i 8 |φ 20 | 2 φ 20 + i 2 (φ 20 -φ 10 ) = 0 (8)
Secular terms in the second of ( 8) are eliminated by requiring that ∂φ 20 ∂τ 0 = 0 ⇒ φ 20 = φ 20 (τ 1 ). Subsequent algebraic manipulation of (8) results in the relation which approximately defines the slow invariant manifold-SIM for 1:1 resonance capture

|φ 10 | 2 =|φ 20 | 2 1 -|φ 20 | 2 2 + λ * 3 4 2 |φ 20 | 4 (9)
This is the manifold wherein the slow dynamics converges in the asymptotic limit τ 0 →∞of the fast dynamics. Hence, the SIM provides us with information regarding the long term behavior of the dynamics under the condition of 1:1 resonance capture. The shape (topology) of the SIM described by ( 9) is depicted in Fig. 9 for the cases of nonlinear damping, λ 3 = 0.3ǫ and absence of nonlinear damping, λ 3 = 0. Numerical simulations were performed on the slow flow equations, system [START_REF] Vakakis | Dynamics of linear discrete systems connected to local essentially nonlinear attachments[END_REF], and compared to direct simulations of the exact system (1) in order to validate the slow flow approximation. The results are depicted in Fig. 10 for two different values of applied impulse I o ; to each simulated response of system (1)we superimpose the corresponding slow variation of the envelope (the slow flow) predicted by system [START_REF] Vakakis | Dynamics of linear discrete systems connected to local essentially nonlinear attachments[END_REF]. In Fig. 10a, the relative response w(t) corresponding to the impulse excitation I o = 0.370 is presented, where the system is locked in 1:1 resonance capture for the full duration of the simulation. The response depicted in Fig. 10c corresponds to a higher impulse value I o = 0.385, where a bifurcation of the dynamics from 1:1 to 1:3 resonance capture takes place. It is clear that the slow flow equations make an excellent prediction for the simulation of Fig. 10a, while the slow flow prediction of the response envelope breaks down in the response of Fig. 10c when the bifurcation from 1:1 to 1:3 resonance capture occurs. This is expected since the slow flow system (6) is only valid under the condition of 1:1 resonance [refer to the ansatz (5)].

In Figs. 10b and10d, the numerical simulations are plotted in comparison with the theoretically predicted SIM for 1:1 resonance capture for the same impulse excitations I o = 0.370 and I o = 0.385, respectively. For the case I o = 0.370, the slow flow simulation follows the lower branch of the SIM for the entirety of the simulation (see Fig. 10b). For I o = 0.385, however, the dynamics initially follows the upper branch of the SIM in the initial phase of the response. That the higher impulse dynamics is initially attracted by the upper branch of the SIM instead of the lower one is due to the fact that in this case the initial energy of the system exceeds the peak energy corresponding to the local fold of the SIM. However, as energy is Fig. 9 Slow invariant manifold (SIM) for 1:1 resonance capture dissipated by (linear and nonlinear) damping the dynamics reaches the local minimum of the SIM, and the response can no longer remain on this SIM. The response subsequently transitions to another slow invariant manifold, which in this particular case is the 1:3 SIM which is discussed later. Since the conditions of 1:1 resonance capture break down after this bifurcation, we need to derive the corresponding slow flow equations under the condition of 1:3 resonance in order to model this transition. This is performed in Sect. 6.

Before we proceed to the study of the bifurcation to 1:3 resonance capture, we note that, following the analysis by Gendelman et al. [START_REF] Gendelman | Dynamics of coupled linear and essentially nonlinear oscillators with substantially different masses[END_REF], we can further simplify the slow flow (6) describing 1:1 resonance capture into a single integral form. This is performed by considering O(ǫ) terms in the slow flow (6)a f t e r applying the multiple scales analysis, and taking the asymptotic limit ǫ→0. Then the O(ǫ) subproblem is derived as

∂φ 10 ∂τ 1 + i + λ * 1 2 φ 10 - i 2 φ 20 = 0 ∂φ 20 ∂τ 1 + λ * 3 -3C * i 8 |φ 20 | 2 φ 20 + λ * 1 2 φ 10 = 0 ( 10
)
We now introduce polar representations for the slow complex amplitudes in order to separate their real moduli and arguments

φ 10 = N 1 exp(iδ 1 ) φ 20 = N 2 exp(iδ 2 ) (11)
which when substituted in ( 10) yield the equivalent slow flow

N ′ 1 e iδ 1 + iN 1 δ ′ 1 e iδ 1 + i + λ * 1 2 N 1 e iδ 1 - i 2 N 2 e iδ 2 = 0 N ′ 2 e iδ 2 + iN 2 δ ′ 2 e iδ 2 ( 12
)
+ λ * 3 -3C * i 8 N 3 2 e iδ 2 + λ * 1 2
N 1 e iδ 1 = 0 where now primes denote differentiation with respect to the slow time scale τ 1 . Dividing the first equation in [START_REF] Lee | Triggering mechanisms of limit cycle oscillations due to aeroelastic instability[END_REF]b ye iδ 1 , dividing the second equation in ( 12) by e iδ 2 , utilizing Euler's representation of a complex number, and setting separately to zero the real and imaginary parts of each equation results in the real slow flow system of equations

N ′ 1 + λ * 1 /2 N 1 -(1/2)N 2 sin δ = 0 N 1 δ ′ 1 + (1/2)N 1 -(1/2)N 2 cos δ = 0 N ′ 2 + λ * 3 /8 N 3 2 + λ * 1 /2 N 1 cos δ = 0 N 2 δ ′ 2 -(3C * /8)N 3 2 + λ * 1 /2 N 1 sin δ = 0 ( 13
)
where δ = δ 1 -δ 2 . These equations are complemented by two additional equations which are derived from the second of ( 8) by imposing the requirement that ∂φ 20 ∂τ 0 = 0, yielding a complex algebraic equation whose real and imaginary parts are set equal to zero, giving

λ * 3 /8 N 3 2 + (N 1 /2) sin δ = 0 (-3C * /8)N 3 2 + N 2 /2 -(1/2) cos δ = 0 (14)
We can now manipulate ( 13) and ( 14) into a single equation for computing the real amplitude N 2

dτ 1 = dη 1 -4η + 3[(λ * 3 /4) 2 + 1]η 2 -(λ * 3 /4)η 2 -λ * 1 {η -2η 2 +[(λ * 3 /4) 2 + 1]η 3 } ( 15
)
where η = N 2 2 . This equation can be integrated by quadratures, so an analytical approximation to the damped transition for 1:1 resonance capture can be derived. Once the slow evolution of the amplitude N 2 (τ 1 ) is analytically approximated from [START_REF] Koumousis | On the dynamic behavior of a lightweight isolator for museum artifacts[END_REF], the other slow varying amplitude N 1 (τ 1 ) and the slow varying phase difference δ(τ 1 ) may be evaluated by means of ( 13) and ( 14). ine the bifurcation to 1:3 resonance capture by again reducing the governing equations of motion by the complexification-averaging technique. However, the new ansatz that will be introduced will necessarily reflect the increased complexity of the dynamics. Based on previous work [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF] capture to 1:3 resonance corresponds to multi-frequency dynamics of the LO and the strongly nonlinear attachment, so resorting to the original system (1) we impose the new ansatz

x ′ + ix = φ 1 (t) exp(it) + φ 3 (t) exp(it/3) v ′ + (i/3)v = φ 2 (t) exp(it/3) + φ 4 (t) exp(it) (16)
It is important to emphasize that the resonant assumption [START_REF] Triplett | Energy harvesting from an impulsive load with essential nonlinearities[END_REF] will not be correct for the general case of sub-harmonic and super-harmonic resonances between the LO and the strongly nonlinear attachment, which would require a different ansatz containing the appropriate harmonic terms. Considering now the expression ( 16), Lee et al. [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF] had shown that in 1:3 resonance capture the response of the LO possesses a dominant frequency component at its linearized natural frequency, i.e., |φ 1 (t)|≫|φ 3 (t)|, whereas the attachment possesses a dominant frequency component at one-third of this frequency, i.e., |φ 2 (t)|≫|φ 4 (t)|. It follows that we can significantly reduce the dimensionality of our study by approximating the response of the system as

x ′ + ix ≈ φ 1 (t) exp(it) v ′ + (i/3)v ≈ φ 2 (t) exp(it/3) (17)
The modified ansatz (17) will be verified by the results of our analysis and will make the study of 1:3 resonance capture analytically possible [indeed, the highdimensionality of ( 16) does not allow a convenient slow flow decomposition of the transient dynamics]. By imposing [START_REF] Quinn | Comparing linear and essentially nonlinear vibration-based energy harvesting[END_REF], we make the assumption that during 1:3 resonance capture the transient response of the LO possesses a single dominant fast frequency equal to unity (its linearized natural frequency), whereas the response of the strongly nonlinear attachment possesses a different dominant fast frequency equal to one third of that.

Substitution of the complex forms (17) into system (1) for the initial conditions (2), and averaging each of the resulting complex equations over their corresponding fast time scales results in the following slow flow that approximately governs the dynamics of 1:3 resonance capture

φ ′ 1 + (λ 1 /2)φ 1 + (1/8)(λ 3 -3Ci)|φ 1 | 2 φ 1 + (9/4)(λ 3 -3Ci)|φ 2 | 2 φ 1 + (9/8)(λ 3 -3Ci)φ 3 2 = 0 ǫφ ′ 2 + (ǫi/6)φ 2 + (9/8)(λ 3 -9Ci)|φ 2 | 2 φ 2 + (1/4)(λ 3 -9Ci)|φ 1 | 2 φ 2 + (3/8)(λ 3 -9Ci)φ 1 φ2 2 = 0 (18)
where overbar denotes complex conjugate. We are now in position to improve our response prediction for the bifurcation from 1:1 to 1:3 resonance capture depicted in Fig. 10c corresponding to I o = 0.385 by directly simulating the slow flow [START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF] in the later part of the transient response. Moreover, since the response of the LO is small compared to the response of the strongly nonlinear attachment in that phase of the response, we will approximately regard the attachment response as being equal to the relative response depicted in Fig. 10c. By extracting the initial conditions from the full simulation at t = 225 and using them to compute the appropriate initial conditions of the slow flow [START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF] at that time instant we obtain the response envelope depicted in Fig. 11. This is in satisfactory agreement with the direct numerical simulation, a fact that validates the previous approximate slow flow analysis for 1:3 resonance capture. Moreover, this improved analytical approximation corrects the result depicted in Fig. 10c. We now focus on a second type of bifurcation that is observed for a lower excitation (energy) level where the system again enters a condition of 1:3 resonance. We will show that in this case 1:3 resonance capture leads to dynamic instability in the response of the nonlinear attachment, manifested as a sudden transition to a relatively high amplitude oscillation. Moreover, it will be interesting to study the response of the attachment prior to the dynamic instability, as it will be shown that it corresponds to a multifrequency decaying oscillation representing a transition between a continuum of different resonance captures.

For I o = 0.277, the time histories of the LO and the strongly nonlinear attachment are presented in Figs. 12a,12b, and the corresponding wavelet spectrum of the relative response between the LO and the 12c. It is remarkable to note that in this case, bifurcation to 1:3 resonance capture is manifested as strong dynamic instability in the response of the dynamics, with the attachment response attaining a relatively high amplitude after the transition to 1:3 resonance capture. An additional remarkable feature of the dynamics is revealed by the wavelet spectrum of Fig. 12d which presents a detail of the wavelet spectrum of Fig. 12c focusing on the initial phase of the dynamics. Indeed, we observe that prior to dynamic instability the dynamics appears to track the manifold of impulsive orbits (IOs) of the underlying Hamiltonian system. This unexpected dynamical behavior implies that the system undergoes multiple resonance captures on multiple resonance manifolds as it tracks the manifold of IOs before finally settling to 1:3 resonance capture and dynamic instability.

This phenomenon is strange and unexpected since, unlike the discrete frequency content of the (countably infinite-Lee et al. [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF]) individual resonance branches (or subharmonic tongues, e.g., S13, S15,... in the Hamiltonian FEP of Fig. 2), the IO manifold possesses orbits with a continuum of frequencies, since it is composed of a countable infinity of resonant periodic IOs and an uncountable infinity of quasiperiodic IOs. Hence, it appears that in the initial phase of the dynamics depicted in Figs. 12a-12d the system undergoes a transition through a continuum of resonance frequencies until the dynamics "locks" into 1:3 resonance capture. This strange phenomenon is robust for lower energy levels as the result depicted in Fig. 12e indicates, corresponding to I o = 0.200. In this case the dynamics is initiated below the energy corresponding to the intersection of the subharmonic tongue S13, so 1:3 resonance capture is not possible. As a result, the dynamics tracks a larger portion of the IO manifold with decreasing energy, and no dynamic instability takes place.

Omitting the initial transition along the IO manifold, whose analytical treatment will be left for a future work, we now analyze the dynamic instability associated with 1:3 resonance capture for the case of impulsive excitation I o = 0.277. In fact, we may divide the region of capture into two distinct phases as depicted in Fig. 13 where the direct numerical simulations are compared to the slow envelopes predicted by the approximate slow flow [START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF]. The slow flow equations were simulated from t = 180 with initial conditions matching those of the direct numerical simulations at that time instant. We note that in the main regime of 1:3 resonance capture the response of the strongly nonlinear attachment attains a nearly constant level of oscillation which is accurately predicted by the slow flow, whereas after the escape from 1:3 resonance capture the envelope of the attachment response decays exactly as predicted by the slow flow. The only region where the slow flow (18) fails to converge with the exact numerical simulation is immediately after the breakup of the dynamic instability (or equivalently in the regime of the dynamics leading to 1:3 resonance capture). This discrepancy is due to the approximations made in deriving the slow flow [START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF], namely the fact that, for the sake of simplicity, we omitted the fast frequency component with normalized frequency (1/3) from the response of the LO and also the fast frequency component with normalized frequency unity from the response of the nonlinear attachment. These simplifications were necessary in order to reduce the dimensionality of the slow flow governing the 1:3 resonance capture, but as Fig. 13 indicates, results in in-accurate modeling of the initial transition period toward dynamic instability.

The near-constant level of amplitude attained by the attachment response in the main regime of 1:3 resonance capture can be analytically predicted by analyzing the stationary response of the undamped slow flow [START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF]. Hence, setting λ 1 = λ 3 = 0 and the time derivatives equal to zero in these equations, we obtain the set of complex algebraic equations

|φ 1 | 2 φ 1 + 18|φ 2 | 2 φ 1 + 9φ 3 2 = 0 - 4ǫ 27C φ 2 + 9|φ 2 | 2 φ 2 + 2|φ 1 | 2 φ 2 + 3φ 1 φ2 2 = 0 (19)
φ 2 =∓|φ 2 | exp( i/3) (20)
where the signs in [START_REF] Manevitch | Complex representation of dynamics of coupled oscillators[END_REF] are the result of multiplying the first equation of ( 19)b y φ1 and recognizing that the first two product terms are real and positive. After substitution of ( 20)into(19), we derive the set of real algebraic equations governing the stationary values of the moduli of φ 1 and φ 2 as

f 1 |φ 1 |, |φ 2 | ≡|φ 1 | 3 + 18|φ 1 ||φ 2 | 2 -9|φ 2 | 3 = 0 f 2 |φ 1 |, |φ 2 | ≡ 2|φ 1 | 2 -3|φ 1 ||φ 2 |+9|φ 2 | 2 - 4ǫ 27C = 0 (21)
It should be clear that the solutions of ( 21) provide the amplitudes of the periodic solutions of the underlying Hamiltonian system (1) corresponding to 1:3 internal resonance between the LO and the attachment as well as the amplitude of the NNM of the Hamiltonian subsystem. In Fig. 14, we depict the solution of (21) as the intersection of two curves that represent the solutions of each of the equations in the (|φ 1 |, 3|φ 2 |) plane. In Fig. 13, this solution is compared to the exact numerical time series of the nonlinear attachment response, and satisfactory agreement is noted.

In addition to the approximated amplitudes of the attachment response at the main phase of 1:3 resonance capture, we also wish to provide an analytical approximation of the slow manifold (analogous to the SIM for 1:1 resonance capture) toward which the response of the damped system transiently evolves. As mentioned above, the main phase of 1:3 resonance capture is characterized by the near constant amplitude of the attachment response. Therefore, one can assume that on the slow manifold of 1:3 resonance capture the time derivative of the second equation of set [START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF]i s approximately equal to zero; this provides us with the approximation for the SIM for 1:3 resonance capture given by

g |φ 1 |, |φ 2 | ≡- 4ǫ 27C φ 2 + 9|φ 2 | 2 φ 2 + 2|φ 1 | 2 φ 2 + 3φ 1 φ2 2 ∼ = 0 (22)
Again, applying simple algebraic manipulations to [START_REF] Lee | Suppression of limit cycle oscillations in the Van der Pol oscillator by means of passive nonlinear energy sinks (NESs)[END_REF] one obtains the relation in terms of amplitudes |φ 1 | and |φ 2 | approximating the SIM for 1:3 resonance capture as

2|φ 1 | 2 -3|φ 1 ||φ 2 |+9|φ 2 | 2 - 4ǫ 27C ∼ = 0 (23)
The algebraic relation (23) represents a projection of the four-dimensional slow flow dynamics onto the plane of amplitudes (|φ 1 |, 3|φ 2 |).I nF i g .15, we depict the approximation of the SIM superimposed on 13).

Finally, we can approximate the escape from 1:3 resonance capture corresponding to decaying response of the strongly nonlinear attachment. This is performed by observing that the response of the LO is nearly zero during this phase. Exploiting this observation, we can greatly simplify the slow flow [START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF] and solve for the simple decay response in closed form. Indeed, assuming that the response of the LO is negligible, φ 1 ∼ = 0, system (18) can be reduced to the equation

ǫφ ′ 2 + ǫi 6 φ 2 + 9 8 (λ 3 -9Ci)|φ 2 | 2 φ 2 ∼ = 0 (24)
Then we take the complex conjugate of relation ( 24) and multiply it by φ 2 ; also, we multiply (24)by φ2 and add it to the previous result to obtain the following first order nonlinear differential equation:

ξ ′ + Dξ 2 = 0 ( 25
)
where ξ = φ 2 φ2 =|φ 2 | 2 and D = 9λ 3 /4ǫ.T h es o l ution to the differential equation ( 25) is expressed in closed form as

ξ(t) = φ 2 (t) 2 = (Dt + P) -1 (26)
where P = P(t o ) is determined from the initial conditions at the start of this phase at time instant t = t o . The comparison of the exact simulation results of system (1)f o rI o = 0.277 with the prediction (26)w a s performed for t o = 1756.9, and satisfactory agreement was observed. This result, combined with the previous analytical derivations, provides a full description of the formation and subsequent decay of the dynamic instability in the system under consideration, and confirms that the governing dynamics of this instability is 1:3 resonance capture. As shown, however, by the numerical simulations of Sect. 4, this is only one of the possible instabilities than can be caused by the geometric nonlinearity. Additional instability scenarios involve different types of resonance captures resulting in multifrequency dynamical transitions in the damped responses. These instabilities can be analyzed applying the analytical methodologies discussed herein.

Concluding remarks

In this work, we documented a rather unexpected result, namely that the addition of geometrically nonlinear viscous damping in a system composed of a linear oscillator-LO coupled to an essentially nonlinear, light attachment, can lead to dynamical instability. This instability appears as a sudden build up of the response of the nonlinear attachment as it engages in resonance capture with the LO, and can be regarded as the transient analog of instabilities that occur in systems with energy input such as the Van der Pol oscillator [START_REF] Lee | Suppression of limit cycle oscillations in the Van der Pol oscillator by means of passive nonlinear energy sinks (NESs)[END_REF] or systems undergoing aeroelastic flutter [START_REF] Lee | Triggering mechanisms of limit cycle oscillations due to aeroelastic instability[END_REF] where limit-cycle oscillations appear at certain ranges of the parameters. Analytical approximations of the dynamics governing instability due to 1:3 resonance capture were determined using a methodology based on complexification and averaging, in spite of the strong stiffness and damping nonlinearities of the system considered. These findings show that (even) purely passive dissipative elements in a structure can lead to dynamic instability due to essential nonlinearities caused by the geometry or the kinematics; hence these results can be of considerable practical significance.

The transient instabilities reported in this work result in strong energy transfer from the LO to the light attachment, which indicates potential for various applications of nonlinear damping, such as in broadband vibration suppression and vibration energy harvesting. An additional interesting dynamical phenomenon reported in this work is the existence of a special type of transient orbit of the strongly nonlinear system during which the dynamics "tracks" the impulsive orbit manifold of the underlying Hamiltonian system. Hence, the system appears to engage in a continuum of resonance captures over certain frequency and energy ranges, or in other words to be in continuous resonance scattering with varying frequency. This unexpected result is the focus of current work by the authors.

Appendix: Derivation of geometrically nonlinear terms

Suppose that we have a spring and damper aligned in parallel horizontally as is shown in Fig. 1. With one end pinned and the other free to translate vertically a distance δ = vx, the stretched length z of the spring can be computed as

z = δ 2 + L 2 (A.1)
The axial force in each spring is computed as

F Si = k(z -L), i = 1, 2( A . 2 )
while the force in each damper is

F Di = λ 2 z ′ ,i = 1, 2( A . 3 )
Computing the time derivative of z z ′ = δδ ′ √ δ 2 + L 2 (A.4)

The total force in each parallel spring and damper combination is

F Ti = F Si + F Di (A.5)
The vertical component of the force in each spring and damper combination is computed as

F Vi = F Ti sin θ (A.6)
where θ is the angle between the horizontal direction and the axis of the rotated spring. From geometry, sin θ is expressed as

sin θ = δ √ δ 2 + L 2 (A.7)
The total vertical force is thus computed as

F Vi = k(δ -L sin θ)+ λ 2 δ ′ sin 2 θ (A.8)

Fig. 3

 3 Fig. 3 Impulsive response of the system with nonlinear damping for I o = 1.0: (a) Relative response v(t) -x(t), (b) wavelet spectrum of relative response, (c) wavelet spectrum on the FEP

Fig. 4

 4 Fig. 4 Impulsive response of the system with nonlinear damping for I o = 0.85: (a) Relative response v(t) -x(t), (b) wavelet spectrum of relative response, (c) wavelet spectrum on the FEP

Fig. 7

 7 Fig. 7 Impulsive response of the system with nonlinear damping for I o = 0.65: (a, b) response and wavelet spectrum of the strongly nonlinear attachment response v(t),(c, d) response and wavelet spectrum of the LO response x(t); (---) natural frequency of the LO

Fig. 8

 8 Fig. 8 Simulated response of system (1) with rescaled parameters for I o = 0.385: (a) Strongly nonlinear attachment time series, (b) LO Time Series, (c) damped transition of the relative displacement superimposed to the Hamiltonian FEP showing bifurcation from 1:1 to 1:3 resonance capture

Fig. 10

 10 Fig. 10 Dynamics of 1:1 resonance capture: (a)T i m es e r i e s and slow flow prediction for I o = 0.370, (b)t r a c k i n gb yt h e slow flow (6)o ft h el o w e rb r a n c ho ft h eS I M( 9)f o rI o =

Fig. 11

 11 Fig. 11 Numerical and slow flow predictions of the bifurcation from 1:1 to 1:3 resonance capture for I o = 0.385

Fig. 12

 12 Fig. 12 Lower energy dynamics: (a, b) Time series of the strongly nonlinear attachment and the LO for I o = 0.277, (c) wavelet spectrum of relative response superimposed to the

Fig. 13 3 φ 1

 1331 Fig.[START_REF] Kerschen | Irreversible passive energy transfer in coupled oscillators with essential nonlinearity[END_REF] Comparison of slow flow[START_REF] Mcfarland | Experimental study of non-linear energy pumping occurring at a single fast frequency[END_REF] and exact numerical simulation for I o = 0.277; the amplitude predicted by the stationary solution,[START_REF] Gendelman | Dynamics of coupled linear and essentially nonlinear oscillators with substantially different masses[END_REF], is also depicted

Fig. 14

 14 Fig. 14 Stationary solution (") of the undamped slow flow (18) corresponding to the periodic solution of the underlying Hamiltonian system (1) under condition of 1:3 internal resonance

Taking the Taylor expansion of sin θ

Assuming small angles

the total vertical force in each spring/damper combination is computed as

Thus, for the two spring/damper pairs, the equation for the vertical force can be generalized to

The validity of the small angle approximation can be preserved by selecting a suitably long distance L in the physical setting.