Roberto De Leo
email: roberto.deleo@ca.infn.it

Sergio Demelio
email: sergio.demelio@ca.infn.it

Numerical analysis of solitons profiles in a composite model for DNA torsion dynamics

We present the results of our numerical analysis of a "composite" model of DNA which generalizes a well-known elementary torsional model of Yakushevich by allowing bases to move independently from the backbone. The model shares with the Yakushevich model many features and results but it represents an improvement from both the conceptual and the phenomenological point of view. It provides a more realistic description of DNA and possibly a justification for the use of models which consider the DNA chain as uniform. It shows that the existence of solitons is a generic feature of the underlying nonlinear dynamics and is to a large extent independent of the detailed modelling of DNA. As opposite to the Yakushevich model, where it is needed to use an unphysical value for the torsion in order to induce the correct velocity of sound, the model we consider supports solitonic solutions, qualitatively and quantitatively very similar to the Yakushevich solitons, in a fully realistic range of all the physical parameters characterizing the DNA.

I. INTRODUCTION

It is now almost thirty years that the role of solitons in fundamental DNA functions is under serious investigation (see [1,2] for a long list of references). In particular a lot of efforts have been made lately [START_REF] Roberts | Nobel Lecture[END_REF][START_REF] Huang | [END_REF][5] to study the phenomenon of the so-called base-flipping, namely the complete opening of a narrow segment of DNA, a phenomenon which is thought to be important for fundamental processes such as replication and transcription.

A successful elementary model for DNA rotational modes, introduced by Yakushevich [6], makes use of a double chain of oscillators where every node of each chain is represented by a disc which interacts with their first neighbors on the same chain and with the disc in front of it on the other chain. The homogeneous version of this model succeds in supporting the existence of topological solitons for a wide range of the physical and geometrical parameters, also allowing analytical solutions for a few cases, but it was shown [6] that adding enough inhomogeneities -i.e. taking FIG. 1: A picture (extracted from [7]) showing the phenomenon of base-flipping, namely the complete opening of a DNA's base pair due to the rotation of (at least) one of them by an angle of π. The four circles on the left enclose respectively a GC base pair, a cytosine base and the two sugars they are respectively attached to on the left chain, while on the right it is shown their center of mass; in the opening process the dihedral angle of the four centers of mass increases from 10 0 to 180 0 . It is clear from the picture that in the process of rotating about the DNA axis both the sugar and base planes also rotate about some other axis; to simplify as much as possible the composite model we will disregard this effect and concetrate only on the rotation about DNA's axis. A n,2

B n,1

B n,2 r C n,1 r C n,2 D n,1 D n,2
FIG. 2: Left. (adapted from [6]) A fragment of our "composite" model of the DNA's double chain: at every node of the chain there are four degrees of freedom, two for each base (the {ϕn,i}) and two for each sugar-phosphate group (the {θn,i}).

Right. Detail of a chain node. Every base {n, i} is allowed to rotate about the C1 atom of the corresponding sugar (see fig. 3), represented by the point Bn,i, only by an angle between [-π/2, π/2] because of the physical contraint represented by the sugar pentagon, so that the effective phase space of the model is (S 1 × [-π/2, π/2]) 2n ; this steric constraint is implemented in the model dynamically through an effective potential. See sec. III for an evaluation of all geometrical constants that appear in the pictures above.

into account the geometrical and dynamical differences between the four possible bases for a realistic DNA segment -causes solitons to lose energy and finally stop their motion after a few nodes on the chain. Our final goal, which will not be reached within this paper, is to determine whether or not realistic segments of DNA support the motion of "rotational" solitons. In the composite model proposed in [1], by separating the degrees of freedom of the bases from the backbone-sugar component, we produced a model where the component supporting the solitons, i.e. the backbone-sugar chain, is perfectly homogeneous and the bases inhomogeneities act just as a perturbation of a homogenous system, leaving hope for a longer life of solitons moving on it. We leave to a future paper the study of the profiles of solitons in inhomogeneous chains and their time evolution.

The numerical results we present in this paper support our guess that the composite model represents an important improvement of the simple torsional models for DNA dynamics in the homogeneous approximation: it supports the existence of solitons and moreover these solitons are very close to the corresponding ones in the Yakushevich model. As a bonus, the increase in the geometrical (and consequentely dynamical) detail turned out to be enough to allow the generation of the solitons corresponding to the Yakushevich ones using for the coupling constants values which are compatible with the physical ones [? ].

As a final remark, we point out that all numerical results and estimates of the geometrical and dynamical parameters included in this paper are an improvement and/or an update of the corrisponding ones published in [1], which is fully superseded by this one.

II. MODEL & EQUATIONS

The composite model for DNA is a homogeneous double chain of coupled double pendula which is a natural generalization of a well known model by Yakushevich where, at every node of each chain, the whole group base-sugarphosphate is represented by a single disc centered at the chain's backbone axis: we simply split the group in two distinct discs, one still centered about the backbone axis and representing the sugar-phosphate group and the second rigidly rotating about a fixed point of the former as shown in fig. 2 (note that, in order to keep the geometry of the model as simple as possible, we consider the chain as plane rather than helicoidal).

Referring to fig. 2, the coordinates of the two discs centers and of the extrema of bases, whose distance we use to determine the intensity of the bonds between the base pairs on the same chain node, are the following:

A n,i (-1) i h/2 , 0 , nδ C n,i (-1) i h/2 + (-1) i+1 (R cos θ n,i + (d bs + r) cos(θ n,i + ϕ n,i )) , (-1) i+1 (R cos θ n,i + (d bs + r) cos(θ n,i + ϕ n,i )) , nδ D n,i (-1) i h/2 + (-1) i+1 (R cos θ n,i + (d bs + 2r) cos(θ n,i + ϕ n,i )), (-1) i+1 (R cos θ n,i + (d bs + 2r) cos(θ n,i + ϕ n,i )), nδ

A c c e p t e d m a n u s c r i p t

where h = 2R + 4r + 2d bs + d eq . The dynamical evolution of our mechanical system is determined by the Lagrangian L = T -V , where T is the kinetic energy and V is the interaction potential. The interactions that are relevant for DNA rotational dynamics are five:

1. The torsion (V t ) between next neighbor sugar-phosphate groups on the same chain, representing the torsional elasticity of the backbone. This force is the result of complex molecular interactions at the backbone level and it is to be considered an "effective" term. Following the principle of keeping the potential expressions as simple as possible until some good reason is found to make them more complicated and keeping in mind that θ can vary on the whole circle, we use for it the "physical pendulum" periodic potential

V t = N -1 n=1 2 i=1 K t [1 -cos(∆θ n,i )] (II.1)
where ∆θ n,i = θ n+1,i -θ n,i . Note that relative angles between next neighbor discs never get big, so that it would be safe also to use its harmonic approximation

V t ≃ N -1 n=1 2 i=1 K t /2[∆θ n,i ] 2 . 2.
The stacking (V s ) between next neighbors bases on the same chain, representing the π -π bonds between the rings that constitute the bases. This interaction is much better understood than the previous one and in particular it is clear that only depends on the relative displacement between next neighbor bases, going rapidly to zero together with the overlapping portion of their surface, e.g. like in a Morse-like potential. However, since also in this case the relative angle of any two bases next to each other keeps small, we simplified its expression by considering it a harmonic bond depending on the "xy" distance between the centers of the bases:

V s = N -1 n=1 2 i=1 1 2 K s d 2 xy (C n+1,i , C n,i ) = N -1 n=1 2 i=1 1 2 K s (d bs + r) 2 [α(cos θ n+1,i -cos θ n,i ) + cos(θ n+1,i + ϕ n+1,i ) -cos(θ n,i + ϕ n,i )] 2 + [α(sin θ n+1,i -sin θ n,i ) + sin(θ n+1,i + ϕ n+1,i ) -sin(θ n,i + ϕ n,i )] 2 (II.2)
where α = R/(d bs + r).

3. The pairing (V p ) between bases on opposite chains, representing the ionic bonds which keep the helices together. This is the best understood force among the ones we are considering and, like in the stacking case, it is known to go rapidly to zero a few Angstrom far from the equilibrium position; in this case though the distance between pairs of bases does get big when a base is flipping and therefore a harmonic approximation result rather unphysical. The interaction can be more realistically modelled by a Morse-like potential [8], nevertheless we will produce profiles also in the harmonic approximation in order to compare our results with those in [6], so we will consider both cases:

V p = N n=1 D 1 -e -µ d(Dn,1,Dn,2 ) 2 ≃ N n=1 1 2 K p (d bs + 2r) 2 [β(cos θ n,1 + cos θ n,2 ) + cos(θ n,1 + ϕ n,1 ) + cos(θ n,2 + ϕ n,2 ) -2 -2β] 2 + [β(sin θ n,1 + sin θ n,2 ) + sin(θ n,1 + ϕ n,1 ) + sin(θ n,2 + ϕ n,2 )] 2 (II.3)
where β = R/(d bs + 2r) and K p = 2Dµ 2 (see eq. (III.1)). Note that, in order to simplify the expression of the elongation from the equilibrium position, we made the "contact" approximation d eq = 0, i.e. we disregarded the inter-bases distance in the equilibrium position; we will show numerically in sec. IV that this does not change significantly the solitons profiles. The same is known to hold also in the Yakushevich model.

4. The helicoidal interactions (V h ) between nucleotides, which are mediated by water filaments (Bernal-Fowler filaments). This is the only ingredient of the composite model which is reminiscent of the helicoidal structure

A c c e p t e d m a n u s c r i p t

of DNA; in particular we will consider those being on opposite helices at half-pitch distance, as they are near enough in three-dimensional space due to the double helical geometry. As the nucleotide moves, the hydrogen bonds in these filaments and those connecting the filaments to the nucleotides are stretched and thus resist differential motions of the two connected nucleotides. We will, for the sake of simplicity and also in view of the small energies involved, only consider filaments between the sugar-phosphate groups, thus only the backbone angles will be involved in these interactions. Recalling that the pitch of the helix corresponds to 10 bases in the B-DNA equilibrium configuration we set

V h = N n=1 2 i=1 K h [1 -cos(θ n+5,i+1 -θ n,i )] (II.4)
where the sum i + 1 is meant modulo 2.

5. The sugar wall (V sw ), representing an "effective" interaction which dynamically restricts the range of the base angles ϕ n,i to some interval [ϕ -, ϕ + ] -ruling out this way the possibility of topologically non trivial configurations for them -to represent the steric constraint represented by the sugars, which prevents the corresponding bases from doing a complete circle about them. We set these angles to φ ± = ±π/2 and verified numerically that the profiles in the composite model converge to the corresponding ones in Yakushevich by narrowing more and more the interval [ϕ -, ϕ + ]. In order not to interfere with the dynamics close to the equilibria positions the potential must be as flat as possible close to the zeroes of the ϕ n,i and must rise rather quickly when the ϕ n,i approach ±π/2. A natural choice, which we implemented in [1], would be to use some high power of the tangent function but divergences cause problems in numerical calculations so we use in this paper a high even power of the sine function

V sw = N n=1 2 i=1 K sw sin 2k (ϕ n,i ) (II.5)
where the coupling constant K sw must be taken big enough to prevent the bases from passing through the ±π/2 barrier but also not so big to interfere too much with the dynamics when the ϕ n,i are closer to the equilibrium position.

As for the kinetic energy, a straightforward calculation shows that for the sugar-phosphate group we have

T t = N n=1 2 i=1 1 2 I t θ2 n,i (II.6)
and for the bases

T s = N n=1 2 i=1 1 2 I s φ2 n,i + 2(1 + α cos ϕ n,i ) φn,i θn,i + (1 + 2α cos ϕ n,i + α 2 ) θ2 n,i (II.7)
Note that by putting the bases angles identically equal to zero the Lagrangian

L = T t + T s -V t -V s -V p -V h -V sw
of the composite model reduces, modulo the helicoidal term, exactly to the Yakushevich homogeneous Lagrangian:

L(θ n,i , 0, θn,i , 0) = N n=1 2 i=1 1 2 I t + (1 + α) 2 I s θ2 n,i - N -1 n=1 2 i=1 K t + (1 + α) 2 (d bs + r) 2 K s [1 -cos(∆θ n,i )] - N n=1 (1 + β) 2 (d bs + 2r) 2 K p {2(1 -cos θ n,1 ) + 2(1 -cos θ n,2 ) -[1 -cos(θ n,1 -θ n,2 )]} (II.8)
Now, once an initial condition for all angles is given, the Lagrangian determines completely the evolution of the state as the "trajectory" q(t) = (θ n,i (t), ϕ n,i (t)) extremizing the action I = q(t) L(q(t), q(t))dt. The initial states we are interested in are those ones which give rise to motions which do not change (much) shape, i.e. which move "rigidly", satisfying therefore the "discrete wave" condition q n+k,i (t) = q n,i (t -kδ/v), where δ is the distance between two consecutive nodes. We look at these solutions as "discrete solitons" able to move on the DNA chain. For a discussion about the analytical properties of the model we refer the reader to the paper [2]. Even in the simplest cases it is impossible to find an explicit solution to the Lagrange equations and therefore the system must be analyzed numerically. To this end we use the fact that in our case h = -δ/v is small enough to claim that

qn,i (t) ≃ q n (t -h) -q n (t) -h = - v δ [q n (t + δ/v) -q n (t)] = - v δ ∆ n,i [q(t)]
so that the kinetic energy becomes

T t = N n=1 2 i=1 I t v 2 2δ 2 [∆θ n,i ] 2 T s = N n=1 2 i=1 I s v 2 2δ 2 [∆ϕ n,i ] 2 + 2(1 + α cos ϕ n,i )∆ϕ n,i ∆θ n,i + (1 + 2α cos ϕ n,i + α 2 )[∆θ n,i ] 2
Having discretized time we are left with just one degree of freedom, the discrete space coordinate n, so that what was previously the system lagrangian L = N n=1 2 i=1 L n,i has become now the action (summed over n) of the discrete Lagrangian L n = 2 i=1 L n,i ; the discrete Lagrange equations, by the minimum action principle, are clearly given simply by dL = (∂ θn,i L, ∂ ϕn,i L, ) = 0 and its solutions, subjected to opportune initial conditions, are exactly the profiles of the solitons that we want to make sure to be present in this generalization of the Yakushevich model.

III. PHYSICAL VALUES OF PARAMETERS & DISPERSION RELATIONS

The evaluation of the main geometrical and dynamical physical quantities entering in the composite model is far from being a trivial matter. There are still very few direct measurements, for example, of the stacking energy between bases within DNA and the torsional backbone force -due to quite complex interactions between the sugar and the phosphate atoms -is just an effective term for which it is hard in principle to define a strenght. Also for geometrical parameters things are not better and often in the literature different values for the same physical quantities are used. In this section we are going to give an estimate for the values of the parameters of the composite model.

A. Kinematical parameters

The kinematical parameters include the masses, the momenta of inertia and the geometrical parameters. We opted for a direct estimate of the geometrical ones from PDB data [START_REF]Pdb repository[END_REF] (note, as already mentioned, that the present results supersede the ones we used in [1]) The masses can be easily calculated on the basis of the chemical structure of DNA. Estimate of momenta of inertia is more involved, especially for the sugar-phosphate group, since the evaluation depends on the choice of the rotation axis and, as pointed out earlier, base-flipping rotations are a complex subject. In our evaluations we made the simplest approximation, namely momenta of inertia of sugar rings have been calculated

A c c e p t e d m a n u s c r i p t

6 with respect to rotations about the main axis of DNA and passing by P atom of the backbone (see fig. 3) and inertia momenta of bases have been calculated with respect to rotations about the same axis by C 1 atom of the sugar ring. The geometrical parameters of interest for the composite model are the longitudinal width of bases l b and of the sugar l s , the distances of the bases from the relative sugars d bs and the distance of a base from the relative dual base d eq . We give our estimates for the masses, moments of inertia and the geometrical parameters for the different bases and their mean values in table I. From those data and using the notations R = l s , r = lb /2, h = l s + d bs + lb + d eq /2, h c = l s + d bs + lb (the bar denotes the mean value), we get the average values for the geometrical parameters appearing in our Lagrangian shown in table II. I: Order of magnitude for the basic geometrical parameters of the DNA. Units of measure are: atomic unit for masses m, 1.67 × 10 -47 Kg • m 2 for the inertia momenta I, Angstrom for l (the longitudinal width of bases and sugar), d bs (the distances sugar-base) and deq (the distance at the equilibrium for the pairs AT and GC). These values have been extracted from the sample "generic" B-DNA PDB data [START_REF]Pdb repository[END_REF], provided by the Glactone Project [START_REF]Glacton project[END_REF], and double checked with the data from [START_REF] Drew | Proc. Natl. Acad. Sci. USA[END_REF], that agrees within 5%. As shown in fig. 3, the lengths of bases and sugar were taken by projecting them on the direction passing through the two phosphate atoms of the chain node to best fit the geometry of the model; the (real) measure for the diameters of the bases A, T, G, C, and the sugar and for d bs are respectively (in Å): 4.6, 4.0, 5.7, 4.0, 3.

B. Dynamical parameters

The determination of the four coupling constants (K p , K t , K s and K h ) characterizing the composite model requires some assessments. Their order of magnitude can be estimated by considering the typical energies of hydrogen bonds (K p ) and the experimental results for the torsional rigidity of the chain (K t and K s ). As we will show later, the geometry of our composite model allows -with values for coupling constants within the experimental ranges -to make predictions about dynamical properties of DNA, as the induced optical frequencies and the phonon velocities, of the same order of magnitude of those experimetally observed (conversely, dispersion relations allow to restrict the selection of the coupling constants to those values which fit better the composite model).

Pairing

The coupling constant K p of the pairing potential (II.3) can be determined by considering the typical energy of hydrogen bonds. The pairing interaction involves two (in the A -T case) or three (in the G -C case) ionic hydrogen bonds. We can model the pairing potential with a Morse function like

V p (d) = D[1 -e -µ(d-deq) ] 2 = 1 2 K p (d -d eq ) 2 + O(d 3 ) , (III.1)
where D is the potential depth, d the distance between the points D n,1 and D n,2 (see fig. are given in [START_REF] Chen | [END_REF] and used in [13]. The values

D = 0.040 eV, µ = 4.45 Å-1
are given in [14] and used in [14][15][16]. Finally, the estimates

D AT = 0.050 eV, D GC = 0.075 eV, µ AT = µ GC = 4 Å-1
are given in [17] and used in [17,18]. The values of coupling constants corresponding to these different values for the parameters appearing in the Morse potential range across an order of magnitude:

3.5 N/m ≤ K p := 2Dµ 2 ≤ 38 N/m . (III.2)
In our numerical investigations we will use a value of K p near to the lower bound given in (III.2); that is, we adopt the value K p = 4 N/m -corresponding to D = 0.030 eV, µ = 2 Å-1 -which leads to an optical excitation threshold of q 4 = 2K p /m b ≃ 32 cm -1 (see eq. (III.12)), so to be in agreement with [19].

Stacking

The determination of the torsion and stacking coupling constants is more challenging and it is based on a smaller amount of experimental data. The main information is the total torsional rigidity of the DNA chain C = Sδ, where δ = 3.4 Å is the base-pair spacing and S is the torsional rigidity. It is known [20,21] that

10 -28 J • m ≤ C ≤ 4 • 10 -28 J • m . (III.3)
This information is used e.g. in [22,23], whose estimate is based on the evaluation of the free energy of superhelical winding; this fixes the range for the total torsional energy to be 180 kJ/mol ≤ S ≤ 720 kJ/mol . (III.4)

In our composite model the total torsional energy of the DNA chain has to be considered as the sum of two parts, the base stacking energy and the torsional energy of the sugar-phosphate backbone. In order to extract the stacking coupling constant we use the further information that π -π stacking bonds amount at the most to 50kJ/mol [24,25].

The stacking potential can be described by a Morse potential whose width σ is of the order 1 Å -so that (see eq. (II.2)) the harmonic approximation is justified -and assuming for the stacking energy the highest value possible, the coupling constant value is K s = 2E s /σ 2 ≃ 16.6N/m. The phonon speed induced by this value of K s is c 3 = δ K s /m b ≃ 3km/s (see eq. (III.12)), which is rather close to the the estimate of 1.8km/s ≤ c 1 ≤ 3.5km/s given in [6].

Torsion and helicoidal couplings

Once picked up the stacking component, our estimate for the torsional coupling constant K t falls in the range 130 kJ/mol ≤ K t ≤ 670 kJ/mol . (III.5)

Assuming (see [26]) that K h ≃ K t /25 = 5 kJ/mol, so that c 4 = 2K t /I s (see equations (III.12) and (III.13)), all of these values for K t induce phonon speeds slightly higher with respect to the estimates cited above, between 5 km/s and 11 km/s. For our numerical investigations, to keep the phonon speed as low as possible, we will set K t = 130 kJ/mol. For numerical calculations it is more convenient to renormalize the Lagrangian with respect to some energy value in order to work with dimensionless quantities. To be coherent with the Yakushevich work we renormalized the Lagrangian with respect to the pairing coefficient l p = K p (d bs + 2r) 2 /2 ≃ 223 kJ/mol; the values of renormalized coefficients for coupling constants g t = K t /l p , g s = K s (d bs + r) 2 /2/l p , g h = K h /l p and g sw = K sw /l p are summarized in table III.

Comparing equation II.8 with the potential used by Yakushevich in [6] it is easy to see that the coupling constant values induced on the Yakushevich model are given by

g = 2g t + 4(1 + α) 2 g s ≃ 30 , K = 2(1 + β) 2 g p ≃ 5.8 (III.6)
so that the normalized Yakushevich coupling constant (the torsional one) corresponding to our choice of parameters is g/K ≃ 5.1. Note that, when using the Morse expression for pairing, the constant in front of the pairing is simply the bond energy D = 0.1 kJ/mol (see eq. (II.3)). We decide nevertheless to divide by the same normalizing factor l p = K p (d bs + 2r) 2 /2 = Dµ 2 so that all dimensionless coupling constants are unchanged except for the pairing one which becomes g ′ p = 1/µ 2 ≃ 0.013.

C. Dispersion relations

The linearized version of the composite Lagrangian is

L lin (θ n,i , ϕ n,i , θn,i , φn,i ) = N n=1 2 i=1 1 2 I t θ2 n,i + 1 2 I s φn,i + (1 + α) θn,i 2 - N -1 n=1 2 i=1 1 2 K t [∆θ n,i ] 2 - N -1 n=1 2 i=1 1 2 K s (d bs + r) 2 [∆ϕ n,i + (1 + α)∆θ n,i ] 2 - N n=1 1 2 K p (d bs + 2r) 2 [ϕ n,1 + ϕ n,2 + (1 + β)(θ n,1 + θ n,2 )] 2 - N n=1 2 i=1 1 2 K h [θ n+5,i+1 -θ n,i ] 2 (III.7)
leading to the following equations of motion:

                   I t θn,i + (1 + α)I s [ φn,i + (1 + α) θn,i ] = -K t θ n,i -K s (d bs + r) 2 (1 + α)[ ϕ n,i + (1 + α) θ n,i ] -K p (d bs + 2r) 2 (1 + β) [ϕ n,1 + ϕ n,2 + (1 + β)(θ n,1 + θ n,2 )] -K h [-θ n+5,i+1 + 2θ n,i -θ n-5,i+1 ] I s [ φn,i + (1 + α) θn,i ] = -K s (d bs + r) 2 [ ϕ n,i + (1 + α) θ n,i ] -K p (d bs + 2r) 2 [ϕ n,1 + ϕ n,2 + (1 + β)(θ n,1 + θ n,2 )]
(III.8)

A c c e p t e d m a n u s c r i p t

where θ n,i = -θ n+1,i + 2θ n,i -θ n-1,i and similarly for ϕ. What we are interested in are the conditions for the existence of wave solutions for the equations above, i.e. in the form θ n,j = Θ ω,k,j e i(kδn+ωt) , ϕ n,j = Φ ω,k,j e i(kδn+ωt) , j = 1, 2

In the approximation α ≃ β, namely considering the bases as pointlike, it is easy to simplify a pair of equations by multiplying the second equation by (1 + α) and substracting it from the first one so that we are left with the pair

I t θn,i = -K t θ n,i -K h [-θ n+5,i+1 + 2θ n,i -θ n-5,i+1 ] , i = 1, 2
After summing and subtracting and then imposing the wave ansatz for the angles on the two equations above we find, by imposing the condition for the existence of non-trivial waves, the first pair of dispersion relations:

ω 2 1 = 4(K t /I t ) sin 2 (kδ/2) + 2(K h /I t )[1 + cos(5kδ)] ω 2 2 = 4(K t /I t ) sin 2 (kδ/2) + 4(K h /I t ) sin 2 (5kδ/2) (III.9)
The second pair of the dispersion relations then can be easily extracted from the remaining pair of equations after the change of coordinates ψ n,i = ϕ n,i + (1 + α)θ n,i , the approximation d bs + 2r ≃ d bs + r (which in turn implies α = β) and finally the fact that I s /(d bs + r) 2 is the mass m b of the base:

m b ψn,1 = -K s ψ n,1 -K p [ψ n,1 + ψ n,2 ] m b ψn,2 = -K s ψ n,2 -K p [ψ n,1 + ψ n,2 ]
Restricting the equations on the wave solutions and summing and subtracting them we find the other pair of dispersion relations

ω 2 3 = 4(K s /m b ) sin 2 (kδ/2) ω 2 4 = 4(K s /m b ) sin 2 (kδ/2) + 2K p /m b (III.10)
Physically, the four dispersion relations correspond to the four oscillation modes of the system in the linear regime. The relations involving ω 1 and ω 2 are associated with torsional oscillations of the backbone. In case of ω 1 there is a threshold for the generation of the excitation originating in the helicoidal interaction, whereas the second torsional mode ω 2 has no threshold and is thus also of acoustical type. The relation involving ω 3 describes relative oscillations of the two bases in the chain with respect to the neighboring bases. As ω 3 (k) → 0 for k → 0 there is no threshold for the generation of these phonon mode excitations. The dispersion relation involving ω 4 describes relative oscillations of two bases in a pair. The threshold for the generation of the excitation is now determined by the pairing interaction. The dispersion relations (III.9,III.10) for values of the physical parameters given in the tables I and III are plotted in Fig. III C; there we plot ω/(2πc), where c is the speed of light (we use the convention, widespread in the literature, of measuring frequencies in 2πc units) versus kδ/2.

The four dispersion relations take a simple form if we consider excitations with wavelength λ much bigger then the intrapair distance, i.e λ ≫ δ; this corresponds to the δ → 0 limit. We have then

ω 2 α -c 2 α k 2 = q 2 α , (III.11)
where c α and q α (α = 1 . . . 4) are, respectively, the velocity of propagation (in the limit k ≫ q α ) and the excitation threshold. They are given by

c 1 = δ (K t -25K h )/I t , q 1 = 2 K h /I t ; c 2 = δ (K t + 25K h )/I t , q 2 = 0; c 3 = δ K s /m b , q 3 = 0; c 4 = δ K s /m b , q 4 = 2K p /m b ; (III.12)
Using the values of the parameters given in the tables I and III we have

c 1 ≃ 0 km/s, q 1 ≃ 4.5 cm -1 ; c 2 ≃ 1 km/s, q 2 = 0 ; c 3 ≃ 3 km/s, q 3 = 0 ; c 4 ≃ 3 km/s, q 4 ≃ 32 cm -1 , (III.13)
where c 1 ≃ 0 comes from the fact that we are taking K t ≃ 25 K h (see table III -this of course just means that c 1 is at least an order of magnitude smaller than the other c i , and therefore negligible). Speeds can be converted to base per seconds by dividing each c i by δ = 3.4 Å; excitation thresholds can be converted in inverse of seconds by multiplying each q i by 2πc, where c is the speed of light.

IV. NUMERICAL ANALYSIS

As we pointed out in sec. II, the profiles of the solitons in the composite model are extremals of the function of 4N variables L(θ n,i , ϕ n,i ), where N is the number of nodes of the chain. The interval of variations of the Lagrangian parameters under our study is such that solitons can have a radius of more than a hundred nodes so a reasonable value of N is of the order of 10 3 ; following Yakushevich [6] we mostly used the value N = 2000.

Clearly the extremals -which by the way in our case turn out to be always maxima -of such a complex function of so many variables can only be obtained numerically; to accomplish that, again following Yakushevich, we choose to use the so-called conjugate gradients method. Ready-to-use implementations of this algorithm are available in the main numerical libraries, in particular in Numerical Recipes (NR, [START_REF]Numerical recipes[END_REF]) and in the GNU Scientific Library (GSL, [START_REF]Gnu scientifica library[END_REF]). As a double check, most of the results presented here have been produced using both implementations; all profiles shown in this paper were produced with GSL.

A. Solitons profiles in the Yakushevich model

As a generalization of the sine-Gordon solitons, also the Yakushevich solitons are "relativistic" and in particular admit a limit speed beyond which no solution exists. We verified numerically that, for a fixed choice of the coupling constants, solitons profiles corresponding to all possible speeds vary by less than 10% from each other, so that it is enough to show the profiles for a single value of the speed. All profiles presented here are relative to static solitons -in other words non-trivial equilibria position of the double chain; in this case the kinetic term is zero and the only parameter left in the renormalized Lagrangian is the coupling constant g of the torsion potential.

As starting point for the extremizing algorithm, following Yakushevich, we used

θ n,i = q i π {1 + tanh [M (2n -N )]} (IV.1)
where (q 1 , q 2 ) is the topological type of the soliton and M a parameter that dictates the steepness of the profile which we use to probe different regions of the phase space.

A first important feature of the model is the disappearance of the discrete solitons when the coupling constant g is smaller than some value g min ≃ 10 (see table V), way before the natural threshold represented by when the solitons become so steep that jump from 0 to 2π in a space shorter than the chain step δ. This behaviour is due to the impossibility of balancing between torsion and pairing when g is too small: indeed both contributions "telescope" but the torsion ones are a priori bounded. A simple way to see how this happen is to give a look to the simplest case, namely the one corresponding to topological numbers (1, 1), when the discrete Lagrange equations can be reduced to a sort of "discrete sine-Gordon" equations, namely g(sin ∆θ i -sin ∆θ i-1 ) = K θ i In each picture the thin line correspond to the physical value for the coupling constant g = 21 and the thick one, to show how the soliton profile changes with g, to g = 150; the continous line is relative to the angles θn,1, the dashed one to the angles θn,2. In the (1, 1) pictures the dashed lines are not visible because the angles profiles for the two angles are identical.

Summing term by term the first n equations we get

g(sin ∆θ n -sin ∆θ 0 ) = K(θ n+1 -θ n -θ 1 + θ 0 )
which, considering that θ 0 = 0 and that for n of the order of N/2 = 1000 also θ 1 ≃ 0, reduces to sin ∆θ n = (K/g)∆θ n . Clearly when g decreases the right term increases and when it becomes bigger than 1 there can be no solution anymore. This feature is particularly relevant for our study because the value g = 5.4 induced by the values of the coupling constants we extracted from literature (see previous section) unfortunately falls in the range where no soliton ariseat least for the basic topological numbers (1, 0), (0, 1) and (1, 1) -in the harmonic approximation for pairing. Since nevertheless 5.4 is rather close to the boundary values for the existence of solutions and the estimate of coupling constants is only qualitative, for the harmonic approximation we increased the value by a factor 4. Note though that this behaviour disappears when a more physical potential (e.g. a Morse one) is used for pairing, leading to the conclusion that the feature above belongs to the harmonic version of the model but does not play any role in realistic models of DNA.

In fig. 6 we show the results for the solitons corresponding to the basic topological numbers (1, 0), (0, 1) and (1, 1) both in the harmonic approximation for pairing (used by Yakushevich in [6]) and with a Morse potential for two values TABLE IV: Energies and diameters for soliton profiles shown in fig. 6; diameters are defined as the number of consecutive nodes between the angles θmin = 1/10 and θmax = 2π -θmin.

(1,0) (0,1) (1,1) g0 7.05 7.05 14.7 TABLE V: The transition values g0 for solitons instability (solitons arise only for g > g0) of the (p, q) solitons.

of the coupling constant g, one corresponding to the physical value and one corresponding to a value about an order of magnitude above in order to enhance the detail in the solitons profile (whose diameter increases monotonically with g). As expected, the profiles corresponding to the Morse potential have a bigger radius, by about a whole order of magnitude, than the corresponding harmonic ones. In all cases the numerical solutions appear to be very robust with respect to the initial configuration, so that we can affirm that in all of them there is a unique extremal.

B. Solitons profiles in the composite model

As a generalization of the sine-Gordon model, also in this composite model there is very little difference between profiles corresponding to different speeds [? ]; in particular, we will consider even in this case only stationary profiles.

As starting point for the extremizing algorithm we used the natural generalization of (IV.1):

θ n,i = q i π {1 + tanh [M (2n -N )]} , φ n,i = 0 .
Since the angles φ n,i are bound to only half of a circle, the corresponding field in the continous approximation cannot describe a topologically non-trivial path and therefore we call these coordinates "non-topological" and there is no integer number associated to them, while the topological numbers q i of the θ angles correspond exactly to the numbers in the simpler model. A further similarity between the two systems is the disappearance of the discrete soliton for too little values of the coupling constants, in this case g t and g s (see fig. 7 -the helicoidal term is at least an order of magnitude smaller than them and it can be disregarded). This is why in all profiles relative to the composite model in the harmonic In the (0, 1) case we compare the profiles to the ones we get by increasing the torsional part to values corresponding to g = 150 in Yakushevich model, once by putting the whole contribute in Vt (thin line) and then putting it wholly in Vs (thinner line). In both cases the profiles diameters get bigger, as expected, but no new qualitative feature appears. In the (1, 1) case we compare the profiles to the ones we get by not making the contact approximation (i.e. keeping deq = 0, thinner line) or by enlarging the helicoidal coupling constant (thin line). Even in this case no new feature appears. approximation (fig. 8) we magnified g s by a factor 4, bringing it to g s = 6.5; when a Morse-like potential is used instead (fig. 9) solitons profiles survive even when both torsion and stacking are turned off, so in that case we keep using the correct value g s = 1.8.

In this section we show the profiles of the solitons for the base topological types (1, 0), (0, 1) and (1, 1) for both the harmonic pairing approximation (fig. 8) and the more physical Morse potential (fig. 9) comparing them with some of the several deformations that can be tested to check their stability.

Harmonic approximation

In fig. 8, first we compare the plot of the (1, 0) soliton with the corresponding soliton in the Yakushevich model; the two graphs are within 10% from each other, testifying that the increase in complexity of the model gives us more features without modifying qualitatively the successful results of the model that generalizes.

Then we compare the profiles of the (0, 1) soliton with those corresponding to a torsional coupling an order of magnitude bigger. Since in this model there are two distinct sources for torsion, we choose to show the profiles corresponding to the extremal cases, namely when all torsion is concentrated in the backbone and when it is instead concentrated in the bases; in this case there are bigger differences between the profiles but we can say that, even after magnifying by an order of magnitude the torsion with respect to the pairing, the picture stays qualitatively the same.

Finally we compare the profiles of the (1, 1) soliton with those obtained by not doing the contact approximation (i.e. by considering the equilibrium position of bases being at distance of 3 Å from each other) and those obtained by magnifying the helicoidal contribution by two orders of magnitude. Even in this case the differences in the profiles turned out negligible.

Summarizing the results, the harmonic approximation turned out to be rather solid even in the composite model but it keeps suffering of the same problem already detected in the Yakushevich model, namely there is no non-trivial discrete soliton if the torsional part is "too weak" and unfortunately the physical values of the coupling constants seem to follow inside this non-existence zone, or at the very least very close to its boundary.

A second concern is that the profiles of the non-topological angles, even the ones generated with higher torsional couplings, in correspondence to the soliton rise from 0 to 2π jump from about -π/2 to π/2 within a single node. This on one side casts some doubt on the existence of a continous conterpart, leading to possible problems for their movement evolution on the chain, and on the other side testifies of a quite violent struggle about the points ±π/2, where two very strong non-linear forces (the pairing and the sugar wall) compete with each other; this behaviour is quite unwelcome firstly because it is unphysical, since after the ionic hydrogen bonds break (and this happens when they get a few Angstrom apart) they do not interact anymore, and secondly because such a violent non-linear interaction most likely destabilizes the system.

Morse potential

In fig. 9, first we compare the plot of the (1, 0) soliton with the corresponding soliton in the Yakushevich model; even with the Morse potential the similarity of the profiles obtained for the two models is striking, keeping within 10% from each other.

Then we compare the plot of the (0, 1) soliton with the ones obtained by decreasing the well width parameter µ, i.e. enlarging the well from 1 Å to 4 Å. As expected we get a narrower profile since in the limit µ → 0 the Morse potential reduces to its harmonic approximation.

Finally we show the (1, 1) soliton profiles and compare them with the ones obtained by neglecting altogether both torsional couplings. The profile is of course narrower but it is still non-trivial, as opposite to the case of the harmonic approximation when even at the physical values for the coupling constants the only profiles we obtain are the constant ones.

Summarizing the numerical analysis relative to the Morse potential, we get the same nice features found in the harmonic approximation but we loose its worst defects. It seems rather clear therefore that a serious investigation of DNA's rotational dynamical properties cannot avoid using a Morse-like potential for pairing, i.e. the pairing coupling must be turned completely off after the distance between bases increases by a few Angstrom.

A c c e p t e d m a n u s c r i p t

(1, 0) 9: Profiles of the solitons corresponding to the smallest non-trivial topological numbers for the composite model with a Morse potential for pairing. In the (1, 0) case we compare the profile with the corresponding Yakushevich one to evidence how close are even in this case the topological component in the composite model to their Yakushevich counterpart. In the (0, 1) case we compare the profiles to the ones we get by decreasing the well width normalized parameter µ(d bs + 2r) from 8.8 to 2 (thinner line). As expected, the profile with a smaller µ gets narrower, due to the fact that the radius of the Morse well got bigger and therefore the pairing interaction looks more like a globally harmonic one. In the (1, 1) case we compare the profiles to the ones we get by disregarding completely both torsional couplings. Differently from what happens in the harmonic approximation, when solitons become trivial for small torsions, the profile has the same qualitative features of the one corresponding to the physical coupling constants.

V. CONCLUSIONS

The numerical analysis we have performed shows the existence of solitonic solutions in the composite DNA model introduced by Gaeta et al. in [1]. The profiles of the topological solitons -in particular, the part relating to the topological degree of freedom -of the composite model are both qualitatively and quantitatively very similar to those of the celebrated Yakushevich model. This means that the most relevant (for DNA transcription) and characterizing feature of the nonlinear DNA dynamics present in the Yakushevich model is preserved by considering geometrically more complex and hence more realistic DNA models.

Moreover, the topological soliton profiles of the composite model seem to change very little when either the physical parameters change in a reasonable range or also the shape of the potential modelling the pairing interaction is modified to a more realistic shape with a sole exception, namely the replacement of the pairing Morse potential with In particular, the shapes of the topological solitons depend very weakly on the interchange of torsional and stacking coupling constant. As a byproduct, this provides justifications for the fact that the Yakushevich model, although based on a strong simplification of the DNA geometry in which torsional and stacking interactions are merged together, works nevertheless quite well in describing solitonic excitations.

The "compositeness" of the composite model becomes relevant -and rather crucial -when it comes, on the one hand, to allow the existence of solitons in correspondence to values of DNA's physical parameters of the correct order of magnitude (with respect to the measures and computational estimates of them available in literature) and, on the other hand, to induce quantities related to small amplitude dynamics, such as transverse phonons speed, which fit experimental observations. In other words, the somewhat more detailed description of DNA dynamics provided by the composite model allows it to be effective -with the same parameters -across regimes, and provide meaningful quantities in both the linear and the fully nonlinear regime.

We expect that the model considered here is the simplest DNA model describing rotational degrees of freedom which, with physically realistic values of the coupling constants and other parameters, allows for the existence of topological solitons and at the same time is also compatible with observed values of bound energies and phonon speeds in DNA. We also expect that solitons may move for considerable distances in this model even in presence of realistic inhomogeneities thanks to the fact that the component that supports the solitons motion, i.e. the sugarphosphate group, is homogeneous and therefore, by separating it from the bases, we can consider the inhomogeneities as a perturbative effect. We leave to a future paper the study of inhomogeneities and time evolution in this model.

2 R

 2 + d bs + 2r 2r + d bs + R 2R + 4r + 2d bs + d eq

1 FIG. 3 :

 13 FIG.3:A node of the DNA chain with a AT pair (left) and a GC pair (right). The measures of the distances are meant as the projections of the actual lenghts over the axis connecting the two phosphates of the same chain node. The distance between the two phosphates is h = 18 Å for both base pairs.

6 × 3

 63 10 3 3.0 × 10 3 4.4 × 10 3 2.3 × 10 3 3.3 × 10 3 2.9 × 10 3 l

  3 and 1.5.) R r d bs deq h hc α β 3.1 1.7 1.0 3.0 18 15 1.1 0.7

1 A 5 FIG. 4 :

 154 FIG. 4: Comparison of Morse potential and its correspondig harmonic approximation

A c c e p t e d m a n u s c r i p t 8 4

 8 . Sugar wallIn order to have a barrier flat enough close to the equilibrium position we choose to use the exponent 2k = 100 in II.5. Numerical tests show that a value of at least K sw = 10 4 kJ/mol is needed in order to keep the bases angles within the [-π/2, π/2] interval.

A c c e p t e d m a n u s c r i p t 10 FIG. 5 :

 105 FIG.5: Graph of the dispersion relations (III.9,III.10) in the first Brillouin zone. We plot ωα/(2πc) (c is the speed of light) as a function of kδ/2. The ω1, ω2, ω3, ω4 are represented respectively by the thick continuous, thin dashed, thin continuous and thick dashed line. Units are cm -1 in the vertical axis and radiants in the horizontal axis.

FIG. 6 :

 6 FIG.6: Profiles of the solitons with the lowest non-trivial topological numbers for the Yakushevich model in the harmonic approximation for pairing (left column) and with a Morse pairing potential (right). In each picture the thin line correspond to the physical value for the coupling constant g = 21 and the thick one, to show how the soliton profile changes with g, to g = 150; the continous line is relative to the angles θn,1, the dashed one to the angles θn,2. In the (1, 1) pictures the dashed lines are not visible because the angles profiles for the two angles are identical.

  q2) g E (MJ/mol) D (nodes) g E (MJ/mol) D (nodes)

FIG. 7 :

 7 FIG.7: Region of the (gt, gs) plane where no solitons arise for the (0, 1) topological type.

A c c e p t e d m a n u s c r i p t 13 (FIG. 8 :

 138 FIG.8: Profiles of solitons corresponding to the smallest non-trivial topological numbers for the composite model in the harmonic approximation for pairing. In the (1, 0) case we compare the profile with the corresponding Yakushevich one to evidence how close the topological component in the composite model are to their Yakushevich counterpart. In the (0, 1) case we compare the profiles to the ones we get by increasing the torsional part to values corresponding to g = 150 in Yakushevich model, once by putting the whole contribute in Vt (thin line) and then putting it wholly in Vs (thinner line). In both cases the profiles diameters get bigger, as expected, but no new qualitative feature appears. In the (1, 1) case we compare the profiles to the ones we get by not making the contact approximation (i.e. keeping deq = 0, thinner line) or by enlarging the helicoidal coupling constant (thin line). Even in this case no new feature appears.

A c c e p t e d m a n u s c r i p t

A c c e p t e d m a n u s c r i p t

  its harmonic approximation, which works fine close to the equilibrium position but fails badly when the bases flip.

TABLE II

 II

: Numerical values (in Å) of the geometrical parameters characterizing the composite model. Here h = 2R + 4r + 2d bs + deq is the width of DNA and hc = h| d 0 =0 is its value in the contact approximation.

TABLE III :

 III Values of the coupling constants (above) and of their normalized counterpart (below).

	Kt	Ks	Kp	K h	Ksw
	130 kJ/mol 16.6 N/m 3.5 N/m 5 kJ/mol 2 • 10 4 kJ/mol
	gt	gs	gp	g h	gsw
	0.58	1.62	1	0.02	100

TABLE VI :

 VI Energies and diameters for soliton profiles shown in fig.8 and 9; diameters are defined as the number of consecutive nodes between the angles θmin = 1/10 and θmax = 2π -θmin.

			Harmonic			Morse
	(q1, q2) gs E (MJ/mol) D (nodes) gs E (MJ/mol) D (nodes)
	(1, 0)		21	54		0.7	136
	(0, 1)	6.5	21	54	1.6	0.7	136
	(1, 1)		63	30		1.1	154

Acknowledgements

We gladly thank G. Gaeta for introducing the problem and both him and M. Cadoni for several enlightening discussions on the subject and for readproofing the manuscript.

A c c e p t e d m a n u s c r i p t


B S Lazarov

J J Thomsen

Using high-frequency vibrations and non-linear inclusions to create metamaterials with adjustable effective properties

Keywords: wave propagation, high-frequency excitation, stiffening, metamaterials, nonlinearity PACS

We investigate how high-frequency (HF) excitation combined with strongly nonlinear elasticity may influence the effective properties for low-frequency wave propagation. The HF effects are demonstrated for linear spring-mass chains with embedded non-linear parts. The investigated mechanical systems can be viewed as a 1D model of materials with non-linear inclusions. The presented analytical and numerical results show that effective material properties can be altered by establishing HF standing waves in the non-linear regions of the chain. In addition, it is demonstrated how true static displacements and forces can be created by using HF excitation with structures having asymmetric displacement-force characteristics.

A c c e p t e d m a n u s c r i p t 1 Introduction

The purpose of this article is to investigate some of the possibilities to change the effective material or structural properties for low frequency (LF) wave propagation by using high-frequency (HF) external excitation combined with strong non-linear elastic material behavior. The ability to create dynamical changes in the materials can be utilized in many problems in structural and mechanical engineering, e.g., creating elasto-mechanical filters, waveguides, ultrasonic devices, opto-mechanical components, sound and vibration isolation, health monitoring. In addition studying further the HF effects of non-linear systems can lead to an improved and safer design of mechanical systems and components.

General HF effects in non-linear systems have been studied analytically, numerically and experimentally in many scientific works, e.g., [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Thomsen | Theories and experiments on the stiffening effect of highfrequency excitation for continuous elastic systems[END_REF][START_REF] Thomsen | Some general effects of strong high-frequency excitation: Stiffening, biasing and smoothening[END_REF][START_REF] Thomsen | Using fast vibrations to quench friction-induced oscillations[END_REF][START_REF] Tcherniak | Slow effects of fast harmonic excitation for elastic structures[END_REF]10,[START_REF] Jensen | Buckling of an elastic beam with added high-frequency excitation[END_REF][START_REF] Fidlin | On asymptotic properties of systems with strong and very strong high-frequency excitation[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF].

HF excitation can change various aspects of the system behavior. This work utilizes the so-called stiffening effect [START_REF] Thomsen | Some general effects of strong high-frequency excitation: Stiffening, biasing and smoothening[END_REF], which results in changes of the effective stiffness for low frequency excitation and the natural frequencies of the system. It can also influence the equilibrium states and their stability. A well known example of the stiffening effect is the upside-down equilibrium stabilization of the Kapitza pendulum on a vibrating support [START_REF] Kapitza | Dynamic stability of a pendulum with an oscillating point of suspension (in Russian)[END_REF], which belongs to the class of parametrically excited systems. The effects of HF parametric excitation are well understood, and have been studied analytically and numerically in many articles [START_REF] Jensen | Buckling of an elastic beam with added high-frequency excitation[END_REF][START_REF] Thomsen | Using fast vibrations to quench friction-induced oscillations[END_REF][START_REF] Thomsen | Some general effects of strong high-frequency excitation: Stiffening, biasing and smoothening[END_REF]. Experimental verification for various physical systems is presented in [START_REF] Jensen | Stiffening effects of highfrequency excitation: experiments for an axially loaded beam[END_REF][START_REF] Yabuno | Experimental investigation of a buckled beam under high-frequency excitation[END_REF][START_REF] Thomsen | Theories and experiments on the stiffening effect of highfrequency excitation for continuous elastic systems[END_REF][START_REF] Thomsen | Effective properties of mechanical systems under high-frequency excitation at multiple frequencies[END_REF]. Stiffening can also be observed in externally excited non-linear systems. Analytical and numerical results for calculating the Thomsen).

A c c e p t e d m a n u s c r i p t

averaged motion of an externally excited system are presented in [START_REF] Thomsen | Using strong non-linearity and high-frequency vibrations to control effective mechanical stiffness[END_REF]. A major part of the references cited above addresses single degree of freedom systems.

HF effects for a general class of multi degree of freedom (MDOF) systems have been studied in [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Thomsen | Some general effects of strong high-frequency excitation: Stiffening, biasing and smoothening[END_REF], and later for continuous systems in [START_REF] Thomsen | Theories and experiments on the stiffening effect of highfrequency excitation for continuous elastic systems[END_REF].

It has been suggested in recent works [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF][START_REF] Thomsen | Using nonlinearity and spatiotemporal property modulation to control effective structural properties: Dynamic rods[END_REF] to use modulation of the structural/material characteristics in space and time, in order to create materials and structures with adjustable effective properties. Such materials are known as dynamic materials [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF]. In [START_REF] Lurie | Low frequency longitudinal vibrations of an elastic bar made of a dynamic material and excited at one end[END_REF][START_REF] Sorokin | Analysis and control of vibrations of honeycomb plates by parametric stiffness modulations[END_REF] HF parametric excitation has been applied to control the effective material properties. In [START_REF] Thomsen | Using nonlinearity and spatiotemporal property modulation to control effective structural properties: Dynamic rods[END_REF] by assuming the existence of a HF standing wave in a continuous non-linear systems, the effective material properties for the LF motion (slow motion) are derived analytically.

Compared to the approaches in [START_REF] Lurie | Low frequency longitudinal vibrations of an elastic bar made of a dynamic material and excited at one end[END_REF][START_REF] Sorokin | Analysis and control of vibrations of honeycomb plates by parametric stiffness modulations[END_REF], where distributed actuation and control are used, the advantage of the approach in [START_REF] Thomsen | Using nonlinearity and spatiotemporal property modulation to control effective structural properties: Dynamic rods[END_REF] is the use of a single excitation source for generating the HF field. However, in contrast to the parametric case, where the stiffness change depends on the frequency of the excitation, for the externally excited non-linear system it does not depend on the frequency of the excitation, and is proportional to the response amplitude and the strength of the non-linearity. If the non-linearities are cubic, the change in the effective material stiffness is predicted to be proportional to the square of the HF response amplitude. The same effect for undamped Fermi-Pasta-Ulam chains has also been predicted analytically in [START_REF] Yu | Nonlinear sinusoidal waves and their superposition in anharmonic lattices[END_REF].

The stiffening can be amplified by increasing the amplitude of the HF standing wave and/or by increasing the strength of the non-linearity. Establishing stable finite amplitude waves in a non-linear medium is not a trivial matter in the general case. Numerical and analytical studies of non-linear chains [START_REF] Rodriguez-Achahc | Generalized instability of periodic traveling waves in aharmonic monatomic chains[END_REF][START_REF] Budinsky | Stability of nonlinear modes and chaotic properties of 1D Fermi-Pasta-Ulam lattices[END_REF] show, that the energy level per particle necessary to start chaotic behavior

A c c e p t e d m a n u s c r i p t

increases with fewer particles in the chain, and thus stable standing waves are easier to establish in shorter chains. This property, combined with the fact that the waves with frequencies above a given threshold cannot propagate along linear chains [START_REF] Brillouin | Wave propagation in periodic structures[END_REF], is used in the present work, in order to create a HF standing wave in a finite non-linear region around the location of the HF excitation. In this way some of the problems related to stability of the HF response can be avoided.

The system considered in the present study is shown in Figure 1. Short nonlinear chains shown with bold lines are embedded in a longer linear one, and the HF excitation is explicitly applied in a non-parametric way in the regions with non-linearities. The HF frequency is within the stop band of the linear chain, and the HF response is localized around the location of the applied excitation.

Various asymptotic methods have been used in the past, to analyze the HF effects in linear and non-linear systems, e.g. [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF][START_REF] Thomsen | Theories and experiments on the stiffening effect of highfrequency excitation for continuous elastic systems[END_REF][START_REF] Thomsen | Some general effects of strong high-frequency excitation: Stiffening, biasing and smoothening[END_REF][START_REF] Tcherniak | Slow effects of fast harmonic excitation for elastic structures[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF]. The Method of Direct Separation of Motion (MDSM) is utilized in this article. The method originates from the work [START_REF] Kapitza | Dynamic stability of a pendulum with an oscillating point of suspension (in Russian)[END_REF]. Later it was generalized by Blekhman [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF], and is utilized in [START_REF] Thomsen | Theories and experiments on the stiffening effect of highfrequency excitation for continuous elastic systems[END_REF][START_REF] Thomsen | Some general effects of strong high-frequency excitation: Stiffening, biasing and smoothening[END_REF][START_REF] Thomsen | Effective properties of mechanical systems under high-frequency excitation at multiple frequencies[END_REF][START_REF] Thomsen | Using nonlinearity and spatiotemporal property modulation to control effective structural properties: Dynamic rods[END_REF][START_REF] Thomsen | Using strong non-linearity and high-frequency vibrations to control effective mechanical stiffness[END_REF][START_REF] Jensen | Buckling of an elastic beam with added high-frequency excitation[END_REF][START_REF] Jensen | Stiffening effects of highfrequency excitation: experiments for an axially loaded beam[END_REF][START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF] . The main assumption in MDSM is that the motion can be represented as a superposition of fast and slow oscillations. The equations for the slow motion are obtained by averaging.

The article is separated into two major parts -theoretical analysis and numerical testing. After introducing the equations of motion for the chain shown in Figure 1, the MDSM is used to obtain the equations governing the slow and fast motions. The obtained analytical results and predictions for the stiffening effect are tested later by numerical simulations.

Mechanical model

The considered system, shown in Figure 1, consists of oscillators with mass m, which are interconnected with springs. The oscillators are labeled with an index j, and their displacements around the static equilibrium position are denoted with ũj . The spring restoring forces are derived from a potential which depends only on the relative distance between the neighbor masses. The equation of motion for the j th mass can be written as:

m d 2 ũj d t2 + c 2 dũ j d t - dũ j-1 d t - dũ j+1 d t + k ∂ Ṽ ∂ ũj = fj (t) (1)
where fj (t) is the external excitation applied to mass j, and c is a linear viscous damping coefficient, which is assumed to be small. By introducing normalized time t = ωt , where ω = k/m, (1) becomes:

d 2 u j dt 2 + c 2 du j dt - du j-1 dt - du j+1 dt + ∂V ∂u j = f j (t) (2)
For a system with N masses, the equations of motion can be written in a matrix form as:

d 2 u dt 2 + C du dt + Φ (u) = f(t) (3)
where u(t) ∈ R N is a vector with an element j equal to the displacement of the j th mass, and Φ (u) is a vector function of the displacements with elements Φ j = ∂V (u)/∂u j . The potential V is assumed to be in the form of power

A c c e p t e d m a n u s c r i p t

series:

V = m k=1 N -1 j=1 a k,j k + 1 r j (u) k+1 + constant (4)
so that:

Φ j = ∂V (u) ∂u j = m k=1 a k,j-1 r j-1 (u) k -a k,j r j (u) k , j = 1 . . . N -1 (5)
where r j (u) = u j+1 -u j , u 0 = 0 and u N = 0. The sum m k=1 a k,j-1 r j-1 (u) k represents the restoring force of the spring interconnecting masses j and j -1.

The external excitation f(t) is assumed to consist of low-frequency and zero mean high frequency components:

f (ωt, Ωt) = f LF (ωt) + f HF (ωt, Ωt) , Ω >> ω (6)
where the high frequency Ω is much larger than the low one ω. In addition, the HF component is assumed to fulfill the following condition:

τ 0 +2π τ 0 f HF (t s , τ ) dτ = 0 (7)
where τ = Ωt and t s = ωt will be called fast and slow times respectively, and t s is assumed independent from τ in the integration.

Approximate equations governing slow and fast motions

The focus is on the LF system response, and on how the HF excitation affects the relevant system properties. Approximate equations governing the slow motion are obtained by using the Method of Direct Separation of Motion (MDSM) [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF]. The solution of (3)-( 7) is assumed to be in the form:

u = z (t s ) + ψ (t s , τ ) (8)

A c c e p t e d m a n u s c r i p t

where z is the slow or the averaged motion of the system and ψ is a fast oscillatory function. The function ψ (t s , τ ) is assumed to be periodic in τ with period 2π, and is required to have a zero fast time average:

ψ (t s , τ ) τ = 1 2π 2π 0 ψ(t s , τ ) dτ = 0 (9)
By averaging (8) and using [START_REF] Tcherniak | Slow effects of fast harmonic excitation for elastic structures[END_REF], the averaged solution u (t s , τ ) τ is equal to the slow motion function z (t s ). The averaging • τ is a linear operator, and given that ( 9) is fulfilled, the following relations are fulfilled as well:

ψ(t s , τ ) ′ τ = ψ(t s , τ ) ′′ τ = ψ(t s , τ ) ′ τ = 0 ( 10
)
where {•} ′ denotes differentiation with respect to τ , and { •} denotes differentiation with respect to t s .

By inserting [START_REF] Thomsen | Using strong non-linearity and high-frequency vibrations to control effective mechanical stiffness[END_REF] and ( 6) into the equations of motion (3), and averaging with respect to the fast time, the equations for the slow motion are obtained in the following form:

ω 2 z + ωC ż + Φ (z + ψ) τ = f LF (t s ) (11)
The expression for the averaged restoring force Φ (z + ψ) τ in ( 11) is obtained by using ( 5), and by expanding the result on LF and HF displacement differences:

Φ j τ = Φ j (z) + Φ j (ψ) τ + m-1 k=1 m p=k+1 c k,p,j-1 r j-1 (z) k + c k,p,j r j (z) k (12) c k,p,j = p k a p,j r j (ψ) p-k τ (13)

A c c e p t e d m a n u s c r i p t

Eq. ( 12) can also be written in the following form:

Φ j τ = Φ eq j (z) + Φ j (ψ) τ (14) Φ eq j (z) = m k=1 a eq k,j-1 r j-1 (z) k + a eq k,j r j (z) k (15)
a eq k,j = a k,j + m p=k+1 a p,j p k r j (ψ) p-k τ (16)
In the general case the term Φ j (ψ) τ = 0 and on average the HF response correspond to additional static forces in the equations governing the slow motion. The HF motion contributes to the effective stiffness and the coefficients in front of all powers of r j (z). An important case can be distinguished when the restoring forces of the springs are odd functions with respect to the relative displacements r j , and the differences r j (ψ) as functions of the fast time are odd functions with respect to a point in the interval (τ, τ + 2π]. In this case the average of the odd powers r j (ψ) 2k+1 τ , k = 0, 1, 2, . . . , is zero, as well as the resulting coefficients in front of the even powers of r i (z). The equivalent static force Φ j (ψ) τ also has a zero value. The slow motion is similar to the response of the system without HF excitation.

An estimate of the HF response ψ is necessary, in order to evaluate the coefficients in the equations governing the slow system motion. The equations governing the HF motion are obtained by inserting ( 8) into (3) and subtracting [START_REF] Chelomei | Dynamic stability upon high-frequency parametric excitation[END_REF], and read as follows:

dψ dt 2 + C dψ dt + Φ (ψ) -Φ (ψ) τ -V(t s ) = f HF (t s , τ ) (17)
where V(t s ) represents the last two sums in [START_REF] Jensen | Buckling of an elastic beam with added high-frequency excitation[END_REF]. The two systems of integrodifferential equations ( 17) and ( 11) describe the motion of the original system completely. However, their solution is not easier than the original equations given by [START_REF] Thomsen | Theories and experiments on the stiffening effect of highfrequency excitation for continuous elastic systems[END_REF]. Further simplifications can be achieved by making assumptions

A c c e p t e d m a n u s c r i p t

for the magnitude of the response. By inserting the following two expressions for the fast motion derivatives:

d dt ψ = ω ψ(t s , τ ) + Ωψ ′ (t s , τ ) (18)
d 2 dt 2 ψ = ω 2 ψ(t s , τ ) + 2Ωω ψ′ (t s , τ ) + Ω 2 ψ ′′ (t s , τ ) (19)
into [START_REF] Lurie | An Introduction to the Mathematical Theory of Dynamic Materials[END_REF], dividing by Ω 2 , and rearranging the terms, the following equation is obtained:

ψ ′′ + 1 Ω Cψ ′ + 1 Ω 2 [Φ (ψ) -Φ (ψ) τ ]+ ω 2 Ω 2 ψ+2 ω Ω ψ′ + ω Ω 2 C ψ- 1 Ω 2 V = 1 Ω 2 f HF (t s , τ ) (20)
The order of the derivatives of the fast motion with respect to the slow and the fast time is assumed to be the same:

O ψ ∼ O (ψ ′ ) ∼ O ψ′ ∼ O ψ ∼ O (ψ) (21)
Waves with frequency ω propagating along the chain have a relatively long wave length compared to a fast wave with frequency Ω. Therefore the differences r j (z) and their powers have relatively small values, and thus the term

V (t s )
/Ω 2 is assumed to be small. For a given response this assumption and ( 21) can be checked a posteriori. Using the above arguments and the initial assumption that Ω ≫ ω, the following simplified system of equations for the fast motion is obtained:

ψ ′′ + + 1 Ω Cψ ′ + 1 Ω 2 [Φ (ψ) -Φ (ψ) τ ] + O ω Ω = 1 Ω 2 f HF (t s , τ ) (22)
In the case when Ω 2 is much larger than the stiffness-to-mass ratio of the springs interconnecting the masses, and there is a lack of resonance effects, the damping and the restoring forces in ( 22) can be neglected, and the so-called

A c c e p t e d m a n u s c r i p t

inertial approximation [START_REF] Blekhman | Vibrational Mechanics -Nonlinear Dynamic Effects, General Approach, Applications[END_REF] for the HF motion is obtained:

ψ ′′ = 1 Ω 2 f HF (t s , τ ) (23)
A discussion on the applicability of the inertial approximation and a comparison with the results obtained by using ( 22) is given in Section 4 and is discussed more generally in [START_REF] Thomsen | Effective properties of mechanical systems under high-frequency excitation at multiple frequencies[END_REF].

Example 1: Linear plus cubic non-linearity

Here the expression for the restoring force for mass j has the following form:

Φ c j (u) = j i=j-1 a 1,i r i (u) + a 3,i r 3 i (u) (24)
The resulting mechanical model corresponds to the Fermi-Pasta-Ulam chain used in computational physics, e.g. [START_REF] Budinsky | Stability of nonlinear modes and chaotic properties of 1D Fermi-Pasta-Ulam lattices[END_REF]. The equivalent restoring force for the equation governing the slow motion is obtained by using ( 14):

Φ c j (u) τ = Φ c j (ψ) τ + j i=j-1 a eq 1,i r i (z) + a eq 3,i r 3 i (z) (25) a eq 1,i = a 1,i + 3a 3,i r 2 i (ψ) τ ( 26
)
a eq 3,i = a 3,i (27)
It is assumed that ( 9) is fulfilled, and hence the coefficients in front of the even powers of r j (z) will be zero. As it appears, for the slow motion, the HF response results in an increase or a decrease of the linear stiffness of the springs between the masses, depending on the sign of a 3,i . The change is proportional to the degree of the non-linearity and the square of the HF amplitude. Similar result for a continuous system are obtained in [START_REF] Thomsen | Using nonlinearity and spatiotemporal property modulation to control effective structural properties: Dynamic rods[END_REF]. The high order non-linear terms in the restoring force are not affected.

A c c e p t e d m a n u s c r i p t 3.2 Example 2: Essentially cubic non-linearity

Here the coefficient a 1,i in ( 24) is zero, and the equivalent coefficient (26) becomes:

a eq 1,i = 3a 3,i r 2 i (ψ) τ (28)
The HF response in this case results in a significant change of the system behavior. Slow waves with small amplitudes will propagate along the HF excited chain in a similar way as in a linear chain with spring stiffness given by (28).

Example 3: Asymmetric non-linearity

Here the restoring force compared to ( 24) is extended with an even order non-linear term and an additional fifth order term.

Φ j (u) = j i=j-1 a 1,i r i (u) + a 2,i r 2 i (u) + a 3,i r 3 i (u) + a 5,i r 5 i (u) (29)
The equivalent coefficients in [START_REF] Fidlin | Nonlinear Oscillations in Mechanical Engineering[END_REF] become as follows:

a eq 1,i = a 1,i + 3a 3,i r 2 i (ψ) τ + 5a 5,i r 4 i (ψ) τ (30) a eq 2,i = a 2,j + 10a 5,j r 3 i (ψ) τ ( 31
)
a eq 3,i = a 3,i + 10a 5,i r 2 i (ψ) τ (32)
The stiffening effect compared to Example 1 is stronger for a HF response with the same amplitude, due to the additional positive term in (31). In addition, due to the even power in the restoring force, the average Ψ (ψ) = 0 and the HF motion will result in an additional equivalent static force for the slow motion.

A c c e p t e d m a n u s c r i p t 3.4 Analytical predictions -summary

The analytical derivations show, that explicitly applied HF external excitation causes changes in the effective stiffness for a LF wave propagation. The stiffness change is proportional to the amplitude of the relative displacements and to the strength of the nonlinearity. In chains with springs having essentially nonlinear behavior, the HF excitation adds effective linear stiffness for the LF motion. The HF effects in nonlinear chains rely on the establishment of a stable HF periodic motion with a large amplitude of the relative displacements.

Numerical simulations and examples

Numerical simulations are performed with a finite spring-mass chain, using standard numerical integration routines from the SUNDIALS library [START_REF] Hindmarsh | SUNDIALS, suite of nonlinear and differential/algebraic equation solvers[END_REF], to test the above analytical predictions. Alteration of the chain stiffness leads to a change of the wave velocity along it. The stiffening effect can be observed by registration of these changes for a system with applied HF excitation. The slow motion has a relatively large wave length, and the contribution of the highorder displacement differences to the restoring force can be neglected. Thus, the wave propagation in the chain with a HF excitation can be compared with the wave propagation in an equivalent linear system with parameters estimated by using [START_REF] Kapitza | Dynamic stability of a pendulum with an oscillating point of suspension (in Russian)[END_REF]. For shorter waves with large amplitude it might be necessary to take into account the high order terms in the equivalent restoring force for the slow motion, however these cases are outside of the scope of this article. The numerical tests are performed for stiffening non-linearity. In the case of a softening, a decrease of the wave velocity can be observed. Depending The simplest model to illustrate the predicted effects, from the numerical point of view, would be a chain with repeatable non-linear springs and masses. Such a system possesses rich dynamic behavior with multiple static and dynamic equilibrium positions. The system response can easily approach chaotic regime.

It is demonstrated in [START_REF] Budinsky | Stability of nonlinear modes and chaotic properties of 1D Fermi-Pasta-Ulam lattices[END_REF][START_REF] Rodriguez-Achahc | Generalized instability of periodic traveling waves in aharmonic monatomic chains[END_REF], that the energy per mass, necessary for obtaining chaotic behavior, is larger for shorter chains, i.e., for shorter chains a larger relative displacement amplitude is necessary, to observe chaotic motion. Therefore, high energy stable oscillatory response can be established by localizing the HF response in short non-linear parts of the chain. Waves can propagate unattenuated along linear chains within specific bands of frequencies called propagation or pass bands, and attenuate within bands of frequencies called stop bands, attenuation zones or band gaps [START_REF] Brillouin | Wave propagation in periodic structures[END_REF]. By embedding the short nonlinear chain in a longer linear one and applying HF excitation with a frequency within the stop bands of the linear chain, the HF response can be localized in the non-linear inclusions.

The investigated system is a chain with 2000 masses. At both ends the displacements are set to be zero. In every cell of twenty springs and twenty masses, one pair of neighboring springs is elastically nonlinear, while the other springs are linear (Figure 2). The HF response is generated by harmonic excitation with amplitude raising slowly from zero to the final value A HF .

f HF = (1 -exp(t/C)) A HF sin (Ωt) (33)
Initially, large damping is introduced in the system, as dashpots between the masses, and the damping coefficients decrease slowly to a prescribed small value with the same rate as the HF amplitude. In this way the transient HF response is damped faster, and the contribution of the damping forces to the dynamic equilibrium, when a steady state is reached, is small. The fast decay of the transients helps the system to reach a steady state oscillatory motion, rather than a complex quasiperiodic or chaotic response. Once a steady state due to the HF excitation is reached, a Gaussian modulated pulse is generated by prescribing the displacement at the left end of the chain:

u 0 (t) = exp (-α (t -t g )) A w cos (ωt -t g ) ( 34
)
where t g = 3500 + 3T , T = 2π/ω, A w is the pulse amplitude, and ω is the slow frequency. The wave velocity for the slow motion can be estimated by finding the displacement maximum u j max and u i max , and by using the following expression:

v w = j -i t j -t i ( 35
)
where t j and t i are the time instances, when mass j and mass i reach the maximum value of the displacements. For the equivalent linear undamped system the group wave velocity can be calculated as a derivative of the dispersion relation ν (ω) with respect to ω, where ν is the phase velocity of a harmonic wave with frequency ω [START_REF] Brillouin | Wave propagation in periodic structures[END_REF].

Inclusions with linear plus cubic non-linear behavior

The stiffening effect is illustrated in Figure 3. Snapshots of the spatial displacements of two chains, with and without HF excitation, are presented. A Gaussian modulated pulse with frequency ω = 0.025 is generated by using (34) in both chains. The pulse in the chain with the HF excitation propagates faster than the one in a chain without it. The small spikes in the system with HF excitation are due to the localized HF response in the non-linear inclusions.

A plot of the wave velocity versus the amplitude of the HF excitation, for a chain with non-linear inclusions with restoring force given by ( 24), is shown in lations.

If the restoring and damping forces in [START_REF] Budinsky | Stability of nonlinear modes and chaotic properties of 1D Fermi-Pasta-Ulam lattices[END_REF] can be considered to be small compared to the inertial force, the inertial approximation [START_REF] Hindmarsh | SUNDIALS, suite of nonlinear and differential/algebraic equation solvers[END_REF] for the HF motion is obtained. An estimate of the wave velocity by using the inertial approximation and numerical simulation with the non-linear system, is shown in Figure 5. As the HF frequency increases, the inertial approximation gives better results. For lower frequencies the estimation starts to deviate from the numerical simulations when the amplitude of the excitation increases.

Increasing the amplitude of the excitation is leading to an increase of the amplitude of the HF response, and thus, the restoring force term in ( 23) is having large contribution to the force balance. Therefore, for highly non-linear systems, the inertial approximation will not lead to a good prediction, even in the case when resonance effects are avoided. For weakly non-linear systems it leads to very good results.

The wave speed in chains with inclusions having different strengths of the non-linearity is shown in Figure 6. As the inertial term is dominant in the equations of motion governing the HF motion, the increase of the excitation amplitude leads to an increase of the amplitude of the HF motion, therefore, stronger non-linearity will lead to a higher equivalent stiffness and a higher wave speed, as shown in Figure 6.

Inclusions with essential non-linearity

HF excitation in combination with strong non-linearity can lead to qualitative changes. If the inclusions in the chain are essentially non-linear (a 1 = 0), it is observed that a wave propagating along the chain will approach a shock wave after passing through a few non-linear inclusions. If a HF excitation is applied, Fig. 6. The left plot shows the group wave velocity for pulse with base frequency ω = 0.025 obtained by numerical simulation for the same system as in Figure 4 with HF excitation frequency Ω = π, and γ = 100, 50, 25, 10 for plots 1,2,3 and 4.

The second plot shows the equivalent stiffness increase of the non-linear spring with respect to the amplitude of the excitation A HF .

the HF motion generates an equivalent linear stiffness, as shown by (28), and the chain becomes transparent for long waves. Comparison of the wave velocity for an equivalent linear system and numerical simulations with the non-linear system, is shown in Figure 7. The amplitude of the generated long wave in this case is A w = 0.1. The system response does not show any unstable behavior in the considered time/space interval, and the results obtained by using the equivalent linear system coincide very well with those obtained by numerical simulations for the non-linear system. The equivalent stiffness in this case can be significantly altered due to the strong non-linearity.

Inclusions with high-order non-linearities

The results presented in Section 3 and derived with the polynomial form of the potential are valid for a much wider class of systems. The equations ( 4) and ( 5) can be viewed as a Taylor expansion of the actual potential or the restoring force. The wave velocity versus the amplitude of the HF excitation for a system with 2 non-linear inclusions on every twenty springs are shown in Figure 8. The restoring forces of the non-linear springs for the three presented cases are given as:

f p = a 1 2n 1 (1 -r(u)) n - 1 (1 + r(u)) n (36)
f e = a 1 2β
1 -e -βr(u) e -βr(u) -1 -e βr(u) e βr(u) (37)

f c = a 1 r(u) + a 1 7β 2 6 r(u) 3 (38)
The values of the parameters β and n are chosen so that the coefficients in front of the third power in the Taylor expansion are equal. The expression for the cubic restoring force (38) is equal to the first two terms of the expansion into

A c c e p t e d m a n u s c r i p t

Taylor series of (37). For small amplitudes of the HF response the equivalent wave speed in the chain is almost the same for all three restoring forces. As the HF amplitude increases, the high order averaged terms start to have a higher contribution to the equivalent coefficients, e.g. ( 30), and the equivalent system becomes stiffer for stronger non-linearity.

Although the results demonstrated so far are for short non-linear embedded chains consisting of two non-linear springs, the analytical predictions in Section 3 are valid for longer non-linear inclusions, as long as the HF motion is periodic or quasiperiodic, and the condition given by ( 9) is fulfilled. The wave velocity versus the amplitude of the HF excitation in a chain with four and five elements in the non-liner inclusion on every twenty springs is shown in Figure 9. The non-linear restoring force of the spring is modeled by using (37) with β = √ 66, and the equivalent linear stiffness is estimated by using ( 38) and (26). The fast time average is estimated by numerical simulations. The HF forces are applied to all masses in the non-linear inclusion with alternating signs, in order to create higher displacement differences.

Generating a static force by means of HF vibrations

According to the analytical predictions in Section 3, a HF response combined with elastic non-linearity can generate an equivalent static force, if the periodic motion is not an odd function with respect to any point in the time interval equivalent static equilibrium position can be obtained by solving the following algebraic equation for z.

Ψ eq (z) + Ψ (ψ) τ = 0 (39)

The possibility to generate an equivalent static force and a new equivalent static equilibrium position, is verified here by using numerical simulations.

The investigated chain, shown in Figure 10, consists of two linear parts with 9 springs and in the middle, 2 non-linear springs shown in bold. The central mass is excited by a harmonic HF force with a frequency in the stop band of the chain, i.e., above the highest linear natural frequency. The stationary HF response is localized around the location of the HF excitation, and the response of the masses far from the excitation is non oscillatory. Therefore, the new equilibrium is a true static equilibrium, and zero mean HF oscillatory motion can be used for generating true static forces and displacements. It should be pointed out that in the case of sub-harmonic resonance, the external HF excitation may generate response within the pass band of the chain, and then the new equilibrium will be static only in the average sense.

The restoring force for the non-linear springs is derived from the widely used in particle simulations Morse potential (e.g. [START_REF] Rodriguez-Achahc | Generalized instability of periodic traveling waves in aharmonic monatomic chains[END_REF]):

V M (r) = a 1 2β 2 1 -e -βr 2 (40)
f M (r) = dV M (r) dr = a 1 β e -βr 1 -e -βr (41)
The strength of the non-linearity is controlled by the parameter β, and the term 1/2β 2 scales (40), so that the slope at r = 0 of the restoring spring force f M (r) is a 1 . The numerical simulations are performed with excitation frequency π and a restoring force for the non-linear springs given by (41). The resulting pre-stressing force in the linear part of the system is compared to the one obtained from the solution of (39) with equivalent parameters evaluated for the Taylor expansion of (41) given as:

f M (r) ≈ a 1 r - 3β 2 r 2 + 7β 2 6 r 3 - 5β 3 8 r 4 + 31β 4 120 r 5 + O r 6 (42)
The time averaged powers of the HF response ψ n τ are estimated from the numerical integration. The results are shown in the first plot in Figure 11.

The amplitude of the HF response is governed mainly by the inertial term for HF excitation. The asymmetry in the force-displacement diagram is having relatively low influence on the amplitude of the HF motion. The pre-stressing force depends mainly on the fast time average Ψ (ψ) τ , and it is decreasing with the parameter β, as the major contribution to the fast time average comes from the even powers of the relative displacements r. This can be seen clearly in the second plot in Figure 11, where the three lines show the results for β = 2.5, 5.0, 10.0. Compared to the original idea [START_REF] Thomsen | Using nonlinearity and spatiotemporal property modulation to control effective structural properties: Dynamic rods[END_REF] for a continuous nonlinear rod, where a single source of the HF excitation is used, here multiple excitation sources are needed, to control the effective material properties. Investigations of the stability properties of the HF response in the nonlinear inclusions, for particular types of nonlinearities, are necessary, in order to utilize the stiffening effect into practical applications. The HF response depends on the HF force amplitude, therefore, resonance effects can also be exploited, in order to minimize the need for a strong HF excitation. The energy required to change the effective material properties depends on the chosen control strategy, and optimal control solutions will be presented in subsequent articles. The presented effects for a 1D chain with non-linear inclusions and HF excitation can be used as a base for developing more complex one, two and three dimensional material

Fig. 1 .

 1 Fig. 1. Spring mass system.

A c c e p t e d m a n u s c r i p t

  on the HF response amplitude and the type of the softening non-linearity, the HF excitation may also result in destabilization of the equilibrium position.

Fig. 2 .

 2 Fig. 2. Part of spring mass system with non-linear inclusions.

Fig. 3 .

 3 Fig. 3. Plots of the displacement response for part of a chain with 2000 masses at time t = 6400.4 with applied HF excitation (solid line) and without HF excitation (dashed line).

Figure 4 .Fig. 4 .

 44 Figure 4. An increase of the velocity, and hence, an increase of the equivalent stiffness is observed, which is in line with the analytical prediction. The velocity obtained by numerical simulations for the non-linear system coincides very well with the one obtained from numerical simulations and analytical calculations for the equivalent linear system. The equivalent parameters are estimated by using (26) with HF response obtained from the numerical simu-

Fig. 5 .

 5 Fig.5. Comparison of the wave speed obtained by using inertial approximation and numerical simulations for the same system as in Figure4, and HF excitation with frequency Ω = 1.5π, 1.8π, 2π shown with solid, dashed and dotted lines respectively.The values obtained by the inertial approximation are shown with point markers.

Fig. 7 .

 7 Fig. 7. The left plot shows the group wave velocity for a pulse with base frequency ω = 0.025 in a chain with two essentially non-linear inclusions on every 20 springs and HF force (amplitude A HF , frequency Ω = π) applied on the mass between the two non-linear springs. The non-linear springs have the following parametersa 1 = 0.0 and a 3 = 0.5γ, where γ = 200. The solid line is obtained by numerical simulations with the non-linear system, the dashed line -by numerical simulations with the equivalent linear system, and the dotted line -by analytical prediction for the undamped system. The second plot shows the equivalent linear stiffness as a function of the excitation amplitude.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Plots of the group wave velocity for a pulse with a base frequency ω = 0.025 in a chain with two non-linear inclusions on every 20 springs and a HF force (amplitude A HF , frequency Ω = π) applied on the mass between two non-linear springs. The restoring force f of the non-linear springs is shown in the second plot as a function of the relative displacements r(u). Solid, dashed and dotted lines correspond to the restoring forces given by equations (36), (37) and (38) respectively, with parameters n = 12 and β = √ 26

A c c e p t e d m a n u s c r i p t f HF Fig. 10 .

 HF10 Fig. 10. Spring mass system used for generating static force by means of HF excitation.

A c c e p t e d m a n u s c r i p t with corresponding restoring force

Fig. 11 . 5 ConclusionsA c c e p t e d m a n u s c r i p t

 115 Fig. 11. Plots of the pre-stressing force in the linear part of the chain, shown in Figure 10, excited by HF harmonic force with amplitude A HF and frequency π. The solid line shows the force obtained from numerical integration, and the dashed line the force estimated by solving (39) with equivalent coefficients estimated by using (42) and β = √ 66. The second plot shows the pre-stressing force for β = 2.5, 5, 10 with solid, dashed and dotted lines respectively.

  models. The possibility to change the stiffness by a HF excitation in parts of the chain can also be used in the design of adjustable band-gap filters for a LF wave propagation.

Acknowledgment

This work was supported by grant 274-05-0498 from the Danish Research Council for Technology and Production Sciences. The authors wish to thank Professor Jakob Soendergaard Jensen for his suggestions and valuable discussions.

A c c e p t e d m a n u s c r i p t