Franz D Fischer

Dieter

T Daxner

F G Rammerstorfer

Stability of Rod-Shaped Nanoparticles Embedded in an Elastic Matrix

Keywords: biological applications, computational mechanics, mechanical behaviour, nanoparticulate materials stability, embedded fibres, load transfer nanoparticulate materials, biological application, stability, embedded fibres, load transfer

In many biological tissues as well as in some technical materials we find nano-sized rod-shaped particles embedded in a relatively soft matrix. Loss of stability of equilibrium, i.e. buckling, is one of the possible failure modes of such materials. In the present paper different kinds of load transfer between matrix and reinforcing particles, which are typical for rod-shaped nanostructures in biological tissues, are considered with respect to stability of equilibrium. Two regimes of matrix stiffnesses leading to different modes of buckling, and a transition regime in between, have been found: soft matrix materials leading to the so-called "flip mode" (also called "tilt mode") and hard matrix materials resulting in "bending mode" buckling. The transition regime is of particular interest for biological tissues.

Numerical and semi-analytical as well as asymptotic concepts are employed leading to results for estimating the critical load intensities both in the form of closed form solutions and diagrams. The analytical solutions are compared with results of finite element analyses. From these comparisons indications are gained for deciding, which of the different analytical approaches should be chosen for a particular nanostructure configuration in terms of the associated buckling modes.

Introduction and Motivation

This paper is motivated by studying the behaviour of mineralized plate-like or rodshaped inclusions in biological tissues, as treated, e.g., by Baohua et al. [START_REF] Baohua | How do slender mineral crystals resist buckling in biological materials? Phil[END_REF], Gao [START_REF] Gao | Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials[END_REF], Jäger and Fratzl [START_REF] Jäger | Mineralized collagen fibrils: A mechanical model with a staggered arrangement of mineral particles[END_REF], as well as Siegmund et al. [START_REF] Siegmund | Failure of mineralized collagen fibrils: Modeling the role of collagen cross-linking[END_REF]. The size of such inclusions is typically in the nano-range. For instance, platelets in the nanostructure of bone are 2 to 4 nm thick and 20 to 100 nm long. Thus, we are talking about nanoparticles.

Load transfer between the organic matrix and the crystallic inclusion as well as stability considerations of such systems are of particular interest. Such situations can be modeled by a stiff elastic rod embedded in a soft elastic matrix leading to an elastically supported beam model. Analytical and numerical treatments of such systems were published already 30 years ago, see, e.g., [START_REF] Madhav | Buckling of finite beams in elastic continuum[END_REF]. With of composites reinforced by collinear short fibres or whiskers the solution of the stability problem of those reinforcements has obtained a revival within the last decade, see, e.g., the papers by Dekret [START_REF] Dekret | Two-dimensional buckling problem for a composite reinforced with a periodic row of collinear short fibers[END_REF], [START_REF] Dekret | Plane problem of stability of composite material armed by periodical tier of short fibers, parallel arrangement[END_REF], [START_REF] Dekret | Near-surface instability of composite materials weakly reinforced by short fibers[END_REF], and Guz and Dekret [START_REF] Guz | On two models in three-dimensional theory of stability of composite materials[END_REF] as well as the literature cited in these papers. It is interesting to note that the problem of magnetic buckling as it occurs due to phase transformations in micromagnetic dots, see [START_REF] Savel'ev | Magnetic and mechanical buckling: Modified Landau theory approach to study phase transitions in micromagnetic disks and compressed rods[END_REF], can be seen as a special case of the more general situations considered here. Furthermore, a similar behaviour is found in [START_REF] Jochum | Proposal for a long-fibre microbuckling scenario during the cure of a thermosetting matrix[END_REF], where, in particular, instability in the form of the "flip mode", i.e., the rigid body rotation (tilting) mode, is presented from experiments.

In Fig. 1 the wavy buckling mode of a thin strip embedded in a soft matrix is shown. Similar experiments are reported in [START_REF] Jochum | Proposal for a long-fibre microbuckling scenario during the cure of a thermosetting matrix[END_REF], where fibre microbuckling during the cure of a fibre reinforced thermosetting polymer is considered.

Strong efforts have been directed to understand the buckling behaviour of fibres both interacting with their bedding and neighbouring fibres by analytical as well as numerical methods.

A similar problem as the one, which was outlined above for composites, is the interaction of slender mineral crystals within biological materials (e.g., mineralized nanoparticles in a protein or collagen matrix) subjected to compressive loading. An excellent overview both on the (staggered) arrangement of those mineral crystals (which can be considered as elastic rods or beams) and the load transfer was given by [START_REF] Jäger | Mineralized collagen fibrils: A mechanical model with a staggered arrangement of mineral particles[END_REF]. Gao [START_REF] Gao | Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials[END_REF] and Baohua et al. [START_REF] Baohua | How do slender mineral crystals resist buckling in biological materials? Phil[END_REF] have investigated the stability behaviour of such nanoparticles by employing the concept of hinged slender plates loaded in the longitudinal direction, as the solution procedure exists since more than 50 years for the stability analysis of rectangular steel plates. The buckling mode shapes documented later in the present paper indicate, however, that a slender hinged plate model may not be the proper mechanical model for thin, elastic plates embedded in a soft matrix. Due the boundary conditions assumed in [START_REF] Gao | Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials[END_REF] and [START_REF] Baohua | How do slender mineral crystals resist buckling in biological materials? Phil[END_REF] the ends are not allowed to move transversally during buckling and, thus, the above mentioned "flip mode" is completely excluded. Since the considered nanoparticles are assumed to be freely embedded in the surrounding matrix and axially loaded, plate-like particles can, within sufficient accuracy, be treated as beams in the context of the presented considerations. This is the reason why in the present paper the considerations are performed for particles of beam shape, having in mind that the results are applicable to platelets, too.

Due to the relevance of the stability analysis of plate-or rod-shaped nanoparticles in biological tissues, represented by non-uniformly axially loaded beams on elastic foundations, the authors are motivated to provide a physically and mathematically sound concept to further the understanding of this particular buckling problem. Although, from the viewpoint of structural mechanics, the problem of buckling of beams on elastic foundation appears to be rather classical, the solutions presented in this paper for the load transfer typical for nanoparticles in soft biological matrices are original.

Formulation of the Mechanical Problem

Mineral crystals embedded in a biological matrix (compare, e.g., the nanostructure of bone) can, like short fibres in composites, be considered as beams which are deformed both in longitudinal direction and by bending. Taking this fact into account the following observations and assumptions must be reflected in a proper physical model:

• The beam is only supported by an elastic bedding. This means that, in contrast to [START_REF] Gao | Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials[END_REF] and [START_REF] Baohua | How do slender mineral crystals resist buckling in biological materials? Phil[END_REF], no essential or geometrical boundary conditions are imposed on the beam. • The bedding may consist of a matrix to which usually averaged elastic properties (Young's modulus E, Poisson's ratio ν) are assigned. Using a continuum approach for modelling the bedding leads to a dependence of the bedding (or foundation) stiffness on the wavelength of the deformation pattern of the beam, in particular, if short-wave buckling is concerned; see, e.g., [START_REF] Vonach | The effects of in-plane core stiffness on the wrinkling behavior of thick sandwiches[END_REF], Vonach and Rammerstorfer [START_REF] Vonach | Wrinkling of thick orthotropic sandwich plates under general loading conditions[END_REF], Vonach and Rammerstorfer [START_REF] Vonach | A general approach to the wrinkling instability of sandwich plates[END_REF], and [START_REF] Tarasovs | Buckling of a coating strip of finite width bonded to elastic half-space[END_REF]. Because of the fact that for not too stiff beddings the buckling pattern of the embedded beams does not show very short wavelengths, the common concept of a constant foundation stiffness k is followed here relating the reaction of the bedding, p, as a transversal force, i.e., acting in the y-direction, per unit length of the beam to the according displacement w by the prominent Winkler relation as

p = -kw. (1)
The unit of the foundation stiffness

[k] is [Force]/[Displacement]/[Length], i.e. [k] = [F ]/[L] 2 .
Here and in the following [". . . "] stands for unit of ". . . ". Thus, [F ] and [L] stands for unit of force and of length, respectively. It is not trivial to define an appropriate value for k. Therefore, estimates of k are usually applied in the literature. In biological structures the elastic properties of the beams and the bedding may differ significantly, see, e.g., [START_REF] Jäger | Mineralized collagen fibrils: A mechanical model with a staggered arrangement of mineral particles[END_REF], [START_REF] Ji | A study of fracture mechanisms in biological nano-composites via the virtual internal bond model[END_REF], or [START_REF] Liu | Poisson ratio can play a crucial role in mechanical properties of biocomposites[END_REF]. Gao [START_REF] Gao | Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials[END_REF] and Baohua et al. [START_REF] Baohua | How do slender mineral crystals resist buckling in biological materials? Phil[END_REF] considered nanostructures in which the stiff nanoplates are arranged in staggered configuration. Such configurations can be found, e.g, in the nanostructure of bone. In such situations and under the assumption that the weakest nanoparticle buckles first and, thus, independently from the others, the following estimate for the bedding stiffness k can be used

k = 2E b b h b , (2)
where E b is the effective Young's modulus of the bedding, i.e. the soft matrix (e.g., collagen), h b is the thickness of the matrix layer, and b is the width of the beam. . For bone typical values of these nanostructural parameters are: For the mineralized nanoparticles (hydroxyapatite) E ≈ 0. 20 to 100 nm. For the matrix (collagen) typical values are E b = 1. × 10 -3 to 7. × 10 -3 µN/nm 2 and h b = 2 to 5 nm, leading to k = 0.4 × 10 -4 to 7. × 10 -3 µN/nm 2 for a width of b = 1.0 nm. (Regarding the choice of b as unit width see Remark 1 below.) For other situations, in which the beam-like nanoparticles are sparsely distributed within the soft matrix (as, e.g., shown in Fig. 1), the determination of the bedding stiffness k is a bit more complicated. It has to be determined depending on the buckling mode, see, e.g., Vonach and Rammerstorfer [START_REF] Vonach | Wrinkling of thick orthotropic sandwich plates under general loading conditions[END_REF] and Fischer and Gamsjäger [6]. • With respect to the load transfer we assume a unit cell of the composite or the biological structure compressed in the axial x-direction. Since the load must be transmitted from the matrix to the fibres, or plates, the assumption of a constant normal force in the longitudinal direction of the beam is generally not justified. Details about the load transfer are reported in the literature, see, e.g., [START_REF] Jäger | Mineralized collagen fibrils: A mechanical model with a staggered arrangement of mineral particles[END_REF][START_REF] Ji | Mechanical properties of nanostructure of biological materials[END_REF]. However, a very sophisticated model developed by Siegmund et al. [START_REF] Siegmund | Failure of mineralized collagen fibrils: Modeling the role of collagen cross-linking[END_REF] shows clearly that also the recently reported concepts, e.g. in [START_REF] Jäger | Mineralized collagen fibrils: A mechanical model with a staggered arrangement of mineral particles[END_REF] and [START_REF] Ji | Mechanical properties of nanostructure of biological materials[END_REF], are too simplified, mainly due to the fact that in the matrix (in case of bones) slip planes occur, and that the fibres are connected to the matrix only in distinct pinning points. Remark 1: It should be mentioned that the width b can be chosen as unit width without loss of generality. In this case the calculated bedding stiffness can be seen as spring stiffness per unit length and unit width. Since the axial load per unit length, t, can be assumed to be constant along the width of the nanoparticle, b, the choice of b = 1.0 leads to the fact that t becomes an axial force per unit length and unit width, i.e., a shear stress loading. This consideration applies to both beam-like and plate-like nanoparticles.

In order to describe the mechanical problem we first fix a coordinate system in the middle of the beam, compare Fig. 2; the longitudinal coordinate x can assume values of -l/2 ≤ x ≤ l/2, with l being the length of the beam. Furthermore, dimensionless coordinates ξ, defined as ξ = x/ [l/2], are introduced, having a range of -1 ≤ ξ ≤ 1. The beam is only supported by a Winkler bedding, so no kinematical boundary conditions appear at ξ ± 1.

A distributed axial force t(ξ) is applied as loading according to Fig. 2. Per definition, this force has a positive sign, if it acts towards the middle of the beam (x = ξ = 0). More specifically, three load cases are investigated:

(1) A constant distributed axial load along the beam, denoted as load case 'CO', compare Fig. 3(a),

t(ξ) = t. (3)
(2) A distributed axial load, which is increasing linearly towards the centre of the beam, with t = 0 at ξ = ±1, and t = t at ξ = 0, We shall refer to this load case as the centre-dominated load case 'CD', compare Fig. 3(c). (3) A distributed axial load which is decreasing linearly towards the centre of the beam, where it vanishes, that is, t = 0 at ξ = 0, and t = ť at ξ = ±1,

t(ξ) = t [1 -|ξ|] . (4
t(ξ) = ť |ξ|. (5)
We shall refer to this load case as the end-dominated load case 'ED', compare Fig. 3(e).

The axial force distribution yields a compressive normal sectional force distribution N (x), or N (ξ), with

N (ξ) = l t max T (ξ) . (6)
The dimensionless function T (ξ) describes the distribution of the normal force N as a function of the position along the beam axis. The maximum axial load t max corresponds to t, t, or ť depending on the load case. For a given axial load distribution t ξ , T (ξ) can be calculated as

T (ξ) = 1 2 1 |ξ| t ξ t max d ξ. ( 7
)
For the three load cases we obtain

(1) load case 'CO':

T (ξ) = 1 2 (1 -|ξ|), compare Fig. 3(b), (2) load case 'CD': T (ξ) = 1 4 [1 -|ξ|] 2 , compare Fig. 3(d), (3) load case 'ED': Ť (ξ) = 1 4 1 -ξ 2 , compare Fig. 3(f).
For the described loading settings we are now looking for equilibrium displacement states w(x), which are non-trivial (w(x) ≡ 0) and fulfill the differential equation for the deflection w(x) of a beam with a non-constant, compressive normal sectional force N (x). In the most general case, this differential equation takes the form

d 2 dx 2 EJ d 2 w dx 2 + d dx N dw dx + kw = 0, (8)
where the product EJ of Young's modulus E and the sectional moment of inertia J defines the bending stiffness of the beam and k is the foundation stiffness.

Equation ( 8) is a linear ordinary differential equation of fourth order with nonconstant coefficients for which no general solution exists. For this reason, approximative solution techniques are applied for finding critical loads.

3. Methods

Finite Element Model

The Finite Element method was used for providing numerical reference solutions over a wide range of geometrical and stiffness parameters. beams were discretized with Euler-Bernoulli beam elements with cubic interpolation functions, thereby implicitly assuming the treatment of slender beams, andin accordance with the analytical models -effectively neglecting shear deformation.

The beams were composed by 1000 beam elements, each. A convergence study revealed that adding more elements did not increase the practically achievable accuracy of the solution, even for those eigenmodes that had the shortest wavelength of all investigated configurations. Results of the finite element simulations will be discussed in Chapter 4.

Energy Method (Ritz Approach)

For the determination of those critical load intensities, which lead to buckling of the beam, first a necessary condition for stable equilibrium has to be defined. According to the Lagrange-Dirichlet theorem, a system containing only conservative forces is stable if the total potential energy Π tot of the system has a strict minimum. This implies also the positive definiteness of the total potential energy, which is an energy function of the structure-load system (see, e.g., [START_REF] Bažant | Stability of Structures[END_REF]).

The elastic beams considered here are treated as discrete elastic systems, their deflection being described by a finite number of n degrees of freedom q 1 , . . . , q n . Including also a load multiplier λ, we obtain Π tot = Π tot (q 1 , . . . , q n ; λ). We now impose a small, geometrical variation on an equilibrium deformation state, which we describe by the variations δq 1 , . . . , δq n of the kinematic degrees of freedom. The corresponding variation δΠ tot of the total potential energy may be expanded into a Taylor series about the equilibrium state,

δΠ tot = δ 1 Π tot + δ 2 Π tot + . . . (9)
with the first and the second variations given by

δ 1 Π tot = n i=1 ∂Π tot ∂q i δq i , (10)
δ 2 Π tot = 1 2 n i=1 n j=1 ∂Π tot ∂q i ∂q j δq i δq j . (11)
The conditions for mechanical equilibrium are δ 1 Π tot = 0 for arbitrary variations δq i , which corresponds to the condition ∂Π tot /∂q i = 0. The condition for stable equilibrium is given by δ 2 Π tot > 0 for arbitrary variations of the degrees of freedom according to the Lagrange-Dirichlet theorem. On the other hand, loss of stability is predicted as soon as the second variation δ 2 Π tot of the total potential energy ceases to be positive definite during a loading process. According to the Trefftz criterion, this requires that the variation of δ 2 Π tot turns to zero, leading to the following conditions for the stability limit:

n i=1 ∂ δ 2 Π tot ∂(δq i ) δq i = 0, or: ∂ δ 2 Π tot ∂(δq i ) = 0 (12)
for all i.

This study is concerned with linear structural systems, for which the second variation of the potential energy coincides with the potential energy itself, δ 2 Π tot = Π tot (as it holds for homogeneous functions of order 2). Consequently, the critical state condition [START_REF] Ji | A study of fracture mechanisms in biological nano-composites via the virtual internal bond model[END_REF] reduces to a condition that takes the form of an equilibrium condition, namely for all degrees of freedom q i .

∂Π tot ∂q i = 0 (13) 7 F
For the present problem the total potential energy Π tot is the sum of the strain energy stored in the beam and in the foundation, Π B and Π F , respectively, as well as the potential Π L of the distributed load:

Π tot = Π B + Π F + Π L . (14)
The axial shortening of the beam due to the compressive normal force will be neglected in the following. Consequently, only the transverse displacement w(x) enters the potential energy Π B .

The individual contributions are detailed as

• the strain energy due to elastic bending of the beam

Π B = EJ 2 l/2 -l/2 d 2 w dx 2 2 dx = EJ 2[l/2] 3 1 -1 d 2 w dξ 2 2 dξ, (15)
• the energy stored in the elastic foundation

Π F = k 2 l/2 -l/2 w 2 dx = kl 4 1 -1 w 2 dξ, (16)
• and the potential of the axial load, which can be found by calculating the product of the local load intensity t and the axial displacement u(x), or u(ξ), of the corresponding load application point, which can, for positive values of x, or ξ, be found as

u = 1 2 x 0 dw dx 2 dx = 1 l ξ 0 dw dξ 2 dξ . . . for x, ξ ≥ 0 ( 17
)
and integrating along the beam axis. Because of the symmetry of the load distribution and the symmetry of the axial displacements, the integral can be calculated over one symmetry half of the beam without loss of information, if the result is multiplied by a factor of two:

Π L = -2 l/2 0 t(x) u(x) dx = -l 1 0 t(ξ) u(ξ) dξ. (18)
Integration by parts, using Eqs. ( 7) and ( 8), gives alternative forms as

Π L = - l/2 0 N (x) dw dx 2 dx = - 2 l 1 0 N (ξ) dw dξ 2 dξ. ( 19
)
The functions w(x), which satisfy the differential equation ( 8) for the deflection of the beam exactly for the given load cases, are not known. We employ the Rayleigh-Ritz method in order to calculate approximative values for the critical load intensities based on the Trefftz criterion [START_REF] Jochum | Proposal for a long-fibre microbuckling scenario during the cure of a thermosetting matrix[END_REF] for linear systems and trial functions w(x) for the deflections which will be detailed in the following sections.

As a first group of Ritz functions, i.e., trial functions used for an approximative description of the deformed configurations w(x) of the buckled beams, polynomials in x are considered. Buckling modes, which are symmetric with respect to the centre of the beam (x = 0), can be described by polynomials of order o S having non-zero coefficients only for even powers of x,

w(x) = oS/2 i=0 a i x 2i . (20)
Buckling modes, which are antisymmetric with respect to the centre of the beam (x = 0), can be defined by polynomials of order o A having nonzero coefficients only for odd powers of x,

w(x) = int(oA/2)+1 i=1 a i x 2i-1 . (21)
The distinction between symmetrical and antimetrical polynomial functions is natural for the representation of the respective buckling modes and, at the same time, makes it easier to solve the respective characteristic equations.

The polynomial Ritz functions presented above do not perform well for foundations with comparatively high stiffness, which lead to short wavelength buckling. Thus, we introduce another set of admissible basis functions w(ξ) for the displacement in y-direction based on the following transcendental and exponential functions:

• The symmetric Ritz function g S g S (ξ) = cos α 2 ξ - 2 α sin α 2 . ( 22
)
The constant term on the right hand side of ( 22) causes the integral 1 -1 g S dξ = 0 to vanish for all values of the shape parameter, α, which reflects bedding force equilibrium in the direction perpendicular to the beam axis.

• The antimetric Ritz function g

A g A (ξ) = sin α 2 ξ . (23)
• The decaying function g DC , which takes the value of one at the centre of the beam, and decays towards the ends

g DC (ξ) = exp - β 2 |ξ| (24)
with β as a further shape parameter. • The decaying function g DE , which assumes the value of one at the ends of the beam and decays towards the centre

g DE (ξ) = exp - β 2 [1 -|ξ|] . (25)
As Ritz functions we use either w(ξ) = q g S (ξ) or w(ξ) = q g A (ξ) directly. These trial functions, which have one degree of freedom only, namely q, will be referred to by labels 'A' and 'S', respectively, in the following. labels are used; e.g., the label 'ADE' is assigned to the single-degree-of-freedom Ritz function wADE (ξ) = q g A (ξ) g DE (ξ), and 'ADC' is assigned to w ADC (ξ) = q g A (ξ) g DC (ξ).

The shape parameters α and β are implicitly defined by Eqs. [START_REF] Vonach | Wrinkling of thick orthotropic sandwich plates under general loading conditions[END_REF][START_REF] Vonach | A general approach to the wrinkling instability of sandwich plates[END_REF] and Eqs. (24, 25), respectively. The parameter α can be interpreted as inverse wavelength parameter of the buckling mode (the full wavelength L w is given by L w = πl/α).

The parameter β describes the intensity, with which the amplitude of the buckling mode decays either toward the ends of the beam, cf. Eqn. (24), or toward its centre, cf. Eqn. (25).

To demonstrate the good approximative nature of the cosinus Ritz function for certain configurations, Fig. 4a shows the symmetric function g S (solid line) plotted against a buckling eigenmode that was predicted by the finite element analysis for a bedding stiffness of k = 10 -5 and a ratio l/h = 100. The case of a constant distributed load 'CO' was considered. It can be seen that the symmetric function predicts the actual mode shape quite well.

Remark 2: It is important to mention that the unit of k, i.e.

[k], is here and in the following to be understood as being consistent so that the product kl has the same unit as the load intensity t, i.e.

[kl] = [t].

Remark 3: It should be noted that the ADE mode with waves near the ends of the beam, see Eqn. (25), is activated for soft bedding only. However, since the rigid body rotation mode, which will be discussed later, is relevant for soft bedding, this ADE mode is not outlined further here.

Increasing the bedding stiffness for the same beam (load case 'CO') to k = 2.8576× 10 -5 gives rise to an antimetric buckling mode, which is depicted in Fig. 4b. The respective Ritz function corresponds to g A . Again, reasonably good agreement is obtained.

Finally, the Ritz function wADE becomes the appropriate choice for approximating the finite element results for a comparatively high bedding stiffness of k = 10 -1 , as is demonstrated for the case of constant axial load 'CO' in Fig. 4c. The dashed lines represent the decaying envelope functions ±g DC . Again, the Ritz function agrees well with the buckling mode predicted by the finite element method.

Since only one degree of freedom, q, is assumed in the Ritz functions, the evaluation of the Trefftz criterion is simplified, yielding for the linear setting

δ δ 2 Π tot = δ(Π tot (q)) = ∂Π tot ∂q δq = 0. ( 26
)
This can be fulfilled for arbitrary variations δq only if ∂Π tot /∂q = 0. If we use the form Π tot = q 2 Πtot (α, β), where Πtot is no longer a function of q, then we come up with the following condition for the critical state: β, appear with respect to which the axial load intensity t * leading to buckling has to be minimized.

∂Π tot ∂q = 2q Πtot = 0, (27)
Inserting the Ritz functions into Eqn. (28) gives the quantity lt max , for which the minimum has to be found. This quantity reads

lt max = EJ l 2 F 1 (α, β) + kl 2 F 2 (α, β) (29)
with the integral terms

F 1 (α, β) = 4 1 0 w 2 dξ 1 0 w 2 T (ξ) dξ, (30)
F 2 (α, β) = 1 4 1 0 w2 dξ 1 0 w 2 T (ξ) dξ. (31
d 2 f(ξ(x)) dx 2
= 4 l 2 f . In order to better understand the relationship between the critical maximum load intensity t max and the shape parameters α and β, respective diagrams are presented in Figs. 5 and6. These diagrams correspond to the load case 'CO', i.e., constant axial load, and the parameter sets that were already discussed in the context of Fig. 4. Figure 5a shows the function t max (α) with β = α/4 for the case of a low bedding stiffness triggering a low wavelength symmetric buckling mode as shown earlier in Fig. 4a. For a slightly higher bedding stiffness, an antimetric buckling eigenform becomes relevant, which is depicted in Fig. 4b with β = α/4. For this mode, the function t max (α) can be found in Fig. 5b.

Finally, for comparatively high bedding stiffnesses, a short-wavelength sinusoidal mode shape with decaying amplitude is predicted according to Fig. 4c. For this case, the maximum intensity t max depends on both parameters α and β. An attempt for the visualization of t max (α, β) is made in Fig. 6.

Looking for the minimum min(lt max ) yields, with the derivatives ∂/∂α, ∂/∂β set to zero, to two highly nonlinear equations in α, β with the solution pairs ᾱ, β being functions of the two system parameters EJ/l 2 , kl 2 . Corresponding solutions can only be found by a numerical procedure, e.g., by applying MAPLE R [http://www.maple.com].

If one sets β = α, 0 < < 1, hundreds of solution pairs ᾱ, β can be found. Therefore, it is more straight-forward to look at the asymptotic behaviour of the system, i.e., large values of α, α 1, for which simple expressions for F 1 and F 2 are found as given in Table 1.

A numerical study has also shown that β can be approximated by a certain fraction of α, e.g., β ≈ α/4 for load case 'CD' and 'ED', so we finally have

F 1 = f 1 α 2 , F 2 = f 2 /α 2 , (32)
f 1 and f 2 being numerical coefficients. Eqs. (29) to (32) allow for a simple estimate of lt * =min(lt max ) and the according value of α, namely α * , as

α * = l 4 k/ [EJ] 4 f 2 /f 1 , (33)
lt * = k [EJ] 2 f 1 f 2 . (34)
As an example we have for the load case 'CD' and β = α/4 S, AS :

4 f 2 /f 1 = 2, f 1 f 2 = 12, (35)
ADC : 4 f 2 /f 1 = 1.84, f 1 f 2 = 4.45. ( 36
)
In the case of soft bedding (relatively small values of kl 4 / [EJ]) the beam looses its stability primarily by a rigid body rotation. The corresponding displacement function is w RB = ϕx for small angles of rotation ϕ. Then the contribution of the strain energy due to bending, Π B ≈ 0, becomes negligible, and lt * can be found directly by setting up an equation for the equilibrium of moments with respect to the centre (x = 0) of the beam. For small positive angles ϕ this equation takes the form

2 l/2 0 t(x) ϕ x dx -2 l/2 0 k ϕ x 2 dx = 0. ( 37
)
Inserting the respective axial load distributions t(x) in Eqn. (37) allows for the calculation of critical load intensities for the three considered load cases as follows:

• for load case 'CO', i.e., constant distributed axial load,

t * RB = 1 3 k l; (38)
• for load case 'CD', i.e., axial load increasing linearly towards the centre of the beam,

t * RB = k l; (39)
• for load case 'ED', i.e., axial load decreasing linearly towards the centre of the beam,

ť * RB = 1 2 k l. ( 40
)
Buckling modes of beams, in which no bending strain energy is activated but only a change in the potential of the loads acting on the beam, have been recently addressed by Grabovsky and Truskinovsky [START_REF] Grabovsky | The flip side of buckling[END_REF] as "flip side of buckling". It should be noted that such modes may interact with bending controlled modes, yielding a transition from the "flip buckling" (tilting) to Euler-type buckling, see also [START_REF] Sidorenko | Reversible switching of hydrogel-actuated nanostructures into complex micropatterns[END_REF].

Apart from the buckling behaviour of the beam on an elastic foundation, we can also extract information about another buckling mechanism from the presented model, namely the symmetric buckling mode of a beam, the centre of which is constrained against translation or rotation. This buckling mode corresponds to the second buckling mode of the beam on the elastic foundation for the limiting case of vanishing foundation stiffness. This situation is treated in the Appendix.

Results

The results obtained from the finite element simulations comprise reference solutions for the critical load intensities t * , t * , and ť * for the constant, the centredominated and the end-dominated axial load cases, respectively. Furthermore, the approximative solutions obtained with the proposed semi-analytical approaches are compared to this reference solutions, and the relative errors, occurring for different Ritz functions, are discussed in detail.

In the first part of this section, the finite element reference solutions are discussed.

For the load case characterized by a constant axial load, i.e., load case 'CO', the critical load intensities t * (k; h/l) are plotted against the foundation stiffness k in Fig. 7.

The solid lines in Figs. Four different curves pertain to thickness-to-length ratios h/l of 0.01, 0.02, 0.05, and 0.10, respectively. Gray areas under a given curve indicate sections of the curves which pertain to symmetric buckling modes. The single gray area in Fig. 7 demonstrates that for the considered system symmetric buckling modes appear only for one limited foundation stiffness interval.

The dashed lines are finite element predictions for the critical load intensities of these beams associated with the second eigenvalues. For very low bedding stiffnesses, these second eigenvalues belong to symmetric eigenmodes such as those discussed in the Appendix. For vanishing bedding stiffnesses, the critical load intensities corresponding to this second eigenvalue asymptotically approach the values derived for the symmetric buckling mode of a beam which is not constrained by any bedding.

Reading Fig. 7 from left to right reveals that for very low bedding stiffness values the buckling mode corresponds more or less to a rigid body rotation ("flip mode") as described above. The analytical solution for this mode is marked by a straight dotted line in Figs. 7 to 9. As the bedding stiffness is increased, the mode changes from an antimetric one to a symmetric one which can be described well by a cosinus function (see Fig. 4a). Increasing the bedding stiffness further causes the buckling mode to become antimetric again. For the considered load case antimetric buckling modes are relevant for bedding stiffnesses above this value.

For increasing bedding stiffness values the wavelength of the buckling mode decreases. Towards the ends of the beam the amplitude of the wavy buckling eigenform decreases asymptotically towards zero. This characteristic behaviour was already discussed in the context of Fig. 4c. For very stiff beddings the critical load intensity approaches asymptotically a value given by Eqs. (34 to 36).

For larger ratios of beam thickness to beam length, i.e., for more compact beams, the critical load intensity curves are shifted towards higher values, taking into account the increased (relative) bending stiffnesses. Of course, the rigid body mode remains unaffected by the bending stiffness, and the asymptotes for beddings of increasing stiffness are all parallel with an inclination of 1/2 in the double-logarithmic Figure 8 shows the critical load intensity ť * for the end-dominated load case 'ED'.

Most of what has already been said for the case of constant axial load applies also for this load case. The big difference is that this load case favors symmetric buckling modes in several different regimes of the bedding stiffness range, as can be seen by the four gray areas in Fig. 8 marking stiffness ranges that cause the relevant buckling mode for h/l = 0.01 to be symmetric. Comparing Fig. 8 to Fig. 7 reveals that the critical load intensity for the load case 'ED' is higher than the one for load case 'CO'. This can easily be attributed to the higher compressive normal force that a given value of t produces in the centre of the beam as compared to the one related to ť.

The third load case with centre-dominated load distribution, i.e. load case 'CD', is associated with critical load intensities t * that are lower than those of the constant and the end-dominated load case. The corresponding critical intensities are plotted over the bedding stiffness in Fig. 9 for four different ratios h/l. Here, the absence of symmetric eigenmodes is most notable. Evidence for the fact that all critical eigenmodes are antimetric regardless of the bedding stiffness can be found in the observation that the dashed lines corresponding to the second eigenvalue do not intersect the solid lines (with the exception of the one which corresponds to the rigid body mode), which represent the lowest eigenvalue.

So far, only the finite element predictions were considered. In the following, the quality of the semi-analytical predictions based on the energy approach with dif- ferent Ritz functions is discussed. First, we take a look at the performance of the polynomial Ritz functions. As a measure for the accuracy of the semi-analytical solutions the relative error with regard to the reference finite element solution is calculated.

For the case of constant axial load and h/l = 0.01, the relative error is plotted over the bedding stiffness k for polynomials of order up to 10 in Fig. 10a. The result for the linear function (n = 1) corresponds to the rigid body mode. This solution is reasonably accurate for bedding stiffnesses up to the value marking the transition to the symmetric buckling mode. The relative error of solutions obtained with transcendental Ritz functions is shown in Fig. 10b.

From Fig. 10 one can conclude that, as long as rather small bedding stiffnesses are concerned, polynomial Ritz functions of sufficiently high order lead to smaller relative errors than the transcendental ones. For situations with higher bedding stiffnesses the transcendental Ritz functions are better. Thus, for all bedding stiffnesses proper Ritz functions exist, based on Eqs. ( 20) to (25), which lead to sufficiently small errors. For the other load cases treated here, similar conclusions can be drawn from the computed results, but are not presented due to space limitations.

So far, the results discussed in Fig. [START_REF] Gao | Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials[END_REF] to Fig.

compare Eqn. (34) in combination with Tab. 1 and Eqn. (32). For rather small values of the bedding stiffness, for which the "flip mode" is relevant (i.e., the bending stiffness EJ is not relevant), the solutions for the critical load intensities can be calculated simply by using Eqs. (38) to (40).

As an example let us consider the nanostructure of bone, for which the following data can be used (b is chosen with b = 1 nm according to Remark 1; cf. Section 2.): For the mineralized platelets, E = 0.1 µN/nm 2 , h = 2 nm, l = 100 nm, hence, h/l = 0.02; for the collagen matrix, E b = 1. × 10 -3 µN/nm 2 , h b = 4 nm, leading to k = 5. × 10 -4 µN/nm 2 . Load case 'CO' is considered. With these data Eqn. (38) leads for the tilt mode t * RB = 1 3 k l = 1.67 × 10 -2 µN/nm 2 . For the bending mode the asymptotic solution used as an estimate can be obtained by using Eqn. (34) in combination with Tab. 1 and Eqn. (32) with the mode shape ADC as t * = 2.57 × 10 -4 µN/nm 2 (an FE analysis has given t * = 2.87 × 10 -4 µN/nm 2 ), which is smaller than t * RB . Thus, bending buckling is the relevant mode. Using Eqn. (41) leads to λ = 1.. Consequently Fig. 7 can directly be used, where one finds for h/l = 0.02 and k = 5.×10 -4 µN/nm 2 : t * ≈ 3×10 -4 µN/nm 2 which is in reasonable agreement with the FE-result, while the estimate resulting from the asymptotic consideration is a bit farer away.

If the values for h and k of the nanostructure of bone used in the above example are changed to become h = 1 nm and k = 1. × 10 -2 µN/nm 2 , i.e. a rather stiff bedding, the estimate for t * using Fig. 7 in conjunction with Eqn. (41) becomes t * = 4.×10 -4 µN/nm 2 . This is exactly that what the finite element analysis renders, too. If, however, the values for E and l in the above bone example are modified in oder to become E = 0.2 µN/nm 2 , l = 20 nm, leading to h/l = 0.1, it becomes obvious that Fig. 7 in conjunction with Eqn. (41) is no longer applicable. Here neither the tilt mode nor the asymptotic estimate are correct; the finite element analysis has lead to a bending buckling mode with a single half wave and a critical load intensity t * = 2.27 × 10 -3 µN/nm 2 .

The above considerations confirm that the proposed estimate procedures work quite well for sufficiently stiff beddings. However, in the transition regime, i.e. between very soft and very stiff bedding, individual numerical or semi-analytical calculations, as they are described in the paper, have to be performed for getting more accurate results. And for the "tilt mode" (or "flip mode") appearing in cases with very soft matrix the simple formulas Eqs. ( 38)-( 40) lead exact values.

Conclusions

In the paper the problem of the stability of equilibrium of nanoparticles in the form of rods or platelets embedded in a relatively soft matrix is treated as buckling of elastically supported beams under axially acting distributed loads, which are introduced to the beam by the surrounding matrix. Computational and semi-analytical as well as asymptotic methods of mechanics are used. Although the appearance of such systems in biological and technical materials is quite widespread, there have been still open questions. This fact is the background motivation for the investiga- tions presented here.

A variety of system parameters and load distributions are investigated in the present paper, and their influence on critical load intensities is presented both in form of formulas and diagrams, which can be used very generally. The quality of the semi-analytically derived results is checked against finite element results. It is demonstrated that, provided the appropriate trial functions are chosen, the Ritz approach leads to reasonably good results.

While for soft beddings polynomial trial functions of sufficient order are appropriate, for stiffer beddings transcendental trial functions should be used.

Obviously two regimes control the stability behaviour of elastic beams or platelets embedded in a relatively soft elastic matrix:

• For the "flip mode", i.e. tilting or pure rotation mode, the critical load intensity is proportional to the bedding stiffness k, i.e lt * ∝ k, compare Eqs. ( 38)-( 40). This mode is particularly relevant for stiffness relations, which can be found in biological tissues (although this mode is not sufficiently perceived in the corresponding literature; see e.g. [START_REF] Baohua | How do slender mineral crystals resist buckling in biological materials? Phil[END_REF], where this mode is not included). hand, [START_REF] Pflüger | Stabilitätsprobleme der Elastostatik[END_REF] does not give a formula for the end-dominated loading case.

The search for a Ritz function, which gives predictions that are at least as good as those in [START_REF] Pflüger | Stabilitätsprobleme der Elastostatik[END_REF] and can be applied to the end-dominated load case, leads to the choice of a tenth order polynomial for the displacement w(ξ),

w(ξ) = w 0 + w 1 ξ 2 + w 2 ξ 4 + w 3 ξ 6 + w 4 ξ 8 + w 5 ξ 10 . (A1)
Solving Πtot = 0 for the respective load intensity t and searching for the minimum t * = min(t(w 0 , w 1 , . . . w 5 )) give the critical load intensities for symmetric buckling under the three considered load cases as follows:

• for load case 'CO' t * S = 62.70

Figure 1 .

 1 Figure 1. Model experiment showing the post-buckling configuration of a foil strip that is glued between two polymer foam blocks. Compression of the block leads to shear stresses along the interface between foil and foam. The resulting compressive normal force distribution leads to buckling.

Figure 2 .Figure 3 .

 23 Figure 2. The elastic beam of length l and bending stiffness EJ in the undeformed configuration on a linear elastic foundation with stiffness k; loading by distributed axial loads t(x) is indicated.

Figure 4 .

 4 Figure 4. Comparison between the buckling eigenmodes predicted by finite element analysis (line with circles) and approximations by Ritz functions for the case of constant axial load 'CO'; a) for a comparatively soft bedding; b) for a little bit stiffer bedding; c) for a stiff bedding.

Figure 5 .

 5 Figure 5. Load intensity tmax as a function of shape parameter α; a) for the symmetric Ritz function (compare Fig. 4a); b) for the antimetric Ritz function (compare Fig. 4b). Load case 'CO' is considered and β = α/4. [F ] is the unit of force, and [L] is the unit of length.

-Figure 7 .

 7 Figure 7. Critical axial load intensity of beams in reference configurations as functions of bedding stiffness; load case 'CO'. Dashed lines represent the second buckling eigenvalue. The dotted line represents the rigid body mode. The gray area refers to symmetric modes. [F ] is the unit of force, and [L] is the unit of length.

8 .

 8 Critical axial load intensity of beams in reference configurations as functions of bedding stiffness; load case 'ED'. Dashed lines represent the second buckling eigenvalue. The dotted line represents the rigid body mode. Gray areas refer to symmetric buckling modes.

Fi gure 9 .

 9 Critical axial load intensity of beams in reference configurations as functions of bedding stiffness; load case 'CD'. Dashed lines represent the second buckling eigenvalue. The dotted line represents the rigid body mode. (No gray area.) diagrams, Figs. 7 to 9.

Figure 10 .

 10 Figure 10. Relative error in the critical axial load intensity t * (load case 'CO'); a) polynomial Ritz functions; b) transcendental Ritz functions. [F ] is the unit of force, and [L] is the unit of length.

( 9 )

 9 refer to beams in reference configuration with E = E ref = 1, b = b ref = 1, and l = l ref = 10 (with consistent units for these quantities).For other values of (E, b, l) these diagrams can also be used for very large bedding stiffnesses (i.e. for the asymptotic case). Here the ordinate values t * , t * , and ť * , respectively, taken from the diagrams according to the corresponding values of k and h/l have to be multiplied by the dimensionless

Figure A1 .

 A1 Figure A1. Sketch showing how the results obtained for the symmetric buckling mode pertain to an equivalent buckling problem.

EJ l 3 , 21 Page

 321 benchmark beams the relative errors of these estimates compared to the finite element solutions, are on the order of 10 -5 . Sketch showing how the results obtained for the symmetric buckling mode pertain to an equivalent buckling problem. 137x99mm (600 x 600 DPI) Model experiment showing the post-buckling configuration of a foil strip that is glued between two polymer foam blocks. Compression of the block leads to shear stresses along the interface between foil and foam. The resulting compressive normal force distribution leads to buckling. 559x247mm (600 x 600 DPI) The elastic beam of length l and bending stiffness EJ in the undeformed configuration on a linear elastic foundation with stiffness k; loading by distributed axial loads t(x) is indicated. 138x81mm (600 x 600 DPI) Normalized intensity of distributed axial load (left column) and resulting normalized normal force distributions (right column) for the three considered load cases. 215x279mm (600 x 600 DPI) Comparison between the buckling eigenmodes predicted by finite element analysis (line with circles) and approximations by Ritz functions for the case of constant axial load (CO); a) for a comparatively soft bedding; 215x279mm (600 x 600 DPI) Load intensity t_max as a function of shape parameter α; a) for symmetric Ritz function 215x279mm (600 x 600 DPI) Load intensity t_max as a function of shape parameter α; b) for the antimetric Ritz function 215x279mm (600 x 600 DPI) Load intensity t_max as a function of shape parameters α and β. 215x279mm (600 x 600 DPI) Critical axial load intensity of beams in reference configuration -load case 'CO' 215x279mm (600 x 600 DPI) Critical axial load intensity of beams in reference configuration -load case 'ED' 215x279mm (600 x 600 DPI) Critical axial load intensity of beams in reference configuration -load case 'CD' 215x279mm (600 x 600 DPI) Relative error in the critical axial load intensity (load case 'CO'); a) polynomial Ritz functions 215x279mm (600 x 600 DPI) Relative error in the critical axial load intensity (load case 'CO'); b) transcendental Ritz functions 215x279mm (600 x 600 DPI)

com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

  Figure 6. Load intensity tmax as a function of shape parameters α and β for the decaying antimetric Ritz function (compare Fig. 4c). [F ] is the unit of force, and [L] is the unit of length.

	October 23, 2009	17:10	Philosophical Magazine		MicroBuckling˙19Oct09
					0.0028
					0.0026
					0.0024
				t max [F] / [L]	0.0018 0.0020 0.0022
					0.0016
					0.0014
					0.0012	0	10	20	β	30	40	50	0	20	40	60 α	80	100 120
			Mode o F S		Load case 'CO' α 2			Load case 'CD' 3α 2	Load case 'ED' 3α 2 /2
			r F 1 A				α 2						3α 2	3α 2 /2
			F 2	ADC (α 2 + 5β 2 )/2 S 16/α 2 P A 16/α 2 e r e ADC 8/(α	(α 2 + 5β 2 ) 48/α 2 48/α 2	(α 2 + 5β 2 ) 16/α 2 16/α 2
									R e v i e w
															O n
															l y

12 http://mc.manuscriptcentral.

2

+ β 2 ) 16/(α 2 + β 2 ) 16/(α 2 + β 2 )

Table

1

. Approximations for functions F1 and F2 for α 1; the notations S, A, and ADC refer to the buckling mode -compare Eqs. (

22

) to (25)

  7 to 9 represent the dependency of the critical load intensities

	October 23, 2009	17:10	Philosophical Magazine	MicroBuckling˙19Oct09
			F o
			r
			P
			e e r
				R e v i e w
				O n
				l y
				14

•

  For the bending modes the critical load intensity is, in an asymptotic sense, proportional to the square root of the bedding stiffness k, i.e. lt * ∝ √ k, see Eqn. (34), or, more specifically, lt * = c√ kEJ, with c being a coefficient depending on loading and the selected Ritz functions.

	October 23, 2009	17:10	Philosophical Magazine	MicroBuckling˙19Oct09
			F o F o
			r r
			P P
		The value of k, where both above mentioned regimes intersect, can be estimated by e e r e e r
		k = As an example, for load case 'CD' Eqn. (42) gives k = 576 EJ/l 4 . In the neigh-EJ l 4 4f 1 f 2 . (42) bourhood of k a transition regime exists, where the simple procedure according to Eqs. (32) to (36) is not directly applicable but the solution pairs (ᾱ, β) can only be found by a numerical study including the minimum search using Eqs. (29) to (31), see Fig. 6. R e v i e w R e v i e w
				O n O n
				l y l y

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

http://mc.manuscriptcentral.com/pm-pml

Acknowledgement

F. D. Fischer expresses his thanks to E. R. Oberaigner for performing the extensive MAPLE R study. T. Daxner and F. G. Rammerstorfer want to thank R. Exler for his assistance in setting up the buckling experiment.

Appendix A. Symmetric Buckling without Foundation

The critical load multipliers corresponding to the second buckling eigenvalue are shown as dashed lines in Figs. 7,8, and 9. Extending those dashed lines towards decreasing foundation stiffnesses shows that the corresponding eigenvalues asymptotically approach values, which represent the buckling eigenvalues of beams subjected to symmetry conditions imposed on their centre point but otherwise freely suspended in space.

These results can be transfered to the mechanically equivalent problem of the buckling of beams, which are fully constrained on one end and free at the other end. In this case the nonuniform distributed load acts towards the constrained end. Due to the symmetry of the problem, the length of the beam in the modified setting is half of the length l of the original beam. Figure A1 illustrates this approach. In Fig. A1(a), a symmetric buckling mode is enforced by constraining the rotation of the beam at its centre. This symmetric buckling mode corresponds to the second buckling mode of the unconstrained problem. Fixing the position of the rotational constraint in space and reducing the foundation stiffness now to zero yield the configuration depicted in Fig. A1(b).

In order to find approximate solutions for the critical load intensities, leading to symmetric buckling of the free beam, the semi-analytical procedure described earlier is applied. The first choice for the Ritz function was the symmetric function g S , compare Eqn. [START_REF] Vonach | Wrinkling of thick orthotropic sandwich plates under general loading conditions[END_REF]. Unfortunately, the predictions obtained with this trial function are worse than the ones provided by the formulas reported in [START_REF] Pflüger | Stabilitätsprobleme der Elastostatik[END_REF]. On the other 20


Fabio ; Moia

Rafael Ferragut
email: rafael.ferragut@polimi.it

Alberto Calloni

Alfredo Dupasquier

Eugenio Carlos

Alberto Macchi

Jian-Feng Somoza

Nie

Vacancy-solute interactions during ageing of a cold worked Mg-RE alloy

Keywords: magnesium alloys, positron annihilation, precipitation Magnesium alloys, Positron annihilation, Precipitation .manuscriptcentral.com/pm-pml

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

Introduction

Precipitation hardening is known to occur in magnesium alloys containing rare-earth elements during artificial ageing of these alloys at elevated temperatures [START_REF] Rokhlin | Magnesium Alloys Containing Rare Earth Metals: Structure and Properties[END_REF]. This phenomenon effectively contributes to the improvement of mechanical properties of this class of magnesium alloys [START_REF] Polmear | [END_REF][START_REF] Polmear | Light Alloys: From Traditional Alloys to Nanocrystals[END_REF]. The precipitation sequence and the morphology of the hardening precipitate phases are relatively well established for commercial alloys that are based on the Mg-Y-RE (rare-earth) system, e.g. WE43 and WE54 alloys [START_REF] Nie | [END_REF][5][6][7][8][9]. An important microstructural factor influencing the nucleation and growth of the strengthening precipitate phases in these age-hardenable alloys is the nature of lattice defects. In the specific case of WE54, which is the alloy selected for study in the present work, it has been reported that cold work after solution treatment and water quench and before artificial ageing at 200°C and at 250°C, can accelerate and enhance the age hardening response [START_REF] Hilditch | Magnesium Alloys and their Applications[END_REF][START_REF] Macchi | [END_REF]. This effect is clearly shown by the hardening curves reported in Ref. [START_REF] Macchi | [END_REF] that are reported in Figure 1 in the present paper for the convenience of readers. The beneficial effect of cold work on age hardening enhancement is also well known for aluminum based alloys (e.g., [12,13]).

Figure 1 about here

In order to gain deeper understanding of the interplay between solute atoms and open volume defects (vacancies and vacancy-like defects that are located in disordered regions and at interfaces) in the WE54 alloy, positron annihilation lifetime spectroscopy has been used by Macchi et al. [START_REF] Macchi | [END_REF]. The propensity of the positrons to be trapped by open volumes inside solids makes positron annihilation lifetime spectroscopy the most sensitive technique for observing vacancy-like defects and the evolution of their density during structural transformation of alloys (for further information on positron annihilation lifetime spectroscopy studies of metallic alloys, see Ref. specific signatures of the chemical species probed by the positrons at the annihilation site [15]. This technique is now widely used for an approximated chemical analysis of vacancy-solute aggregates in light alloys (an example with full details on the methodology is given in Ref. [16]).

The CDB measurements of artificially aged, undeformed samples of the same alloy have been reported in Ref. [17]. These measurements reveal the almost total absence of frozen-in vacancies in the as-quenched condition and suggest that only thermal equilibrium vacancies contribute to solute transport. In the present work, the positron annihilation study is extended to deformed (cold worked) specimens. These additional data and their comparison with the results reported in Refs. [START_REF] Macchi | [END_REF]17] give us the opportunity to shed light on the mechanism in which deformation-induced defects influence precipitation.

Experimental procedures

Sample preparation

The samples were taken from a cast plate of WE54 (Magnesium Elektron Company, England), with a nominal composition of Mg -5.1Y-3.3RE -0.4Zr (wt.%), where RE stands for a mixture of Nd and heavy rare-earth elements, with Nd as the main component. The samples were solution treated for 8 hours at 525°C, quenched in water at room temperature and immediately cold rolled to a thickness reduction of 6% and 12%. After a dwell time of a few days at room temperature (this alloy is not prone to natural ageing at room temperature), the samples were artificially aged at 250°C for up to 4 hours. The ageing treatment was temporarily interrupted during positron annihilation measurements. The interruption entailed water quenching followed by a short exposure at room temperature before 24 hours dwelling at liquid nitrogen temperature. In the text below, the effect of artificial ageing interruption is neglected and the accumulated dwell time at 250°C is taken as the effective ageing time.

Positron annihilation lifetime spectroscopy

Positron annihilation lifetime measurements published in Ref. [START_REF] Macchi | [END_REF] were performed at room temperature. Since the positron trapping rate in the diffusion-controlled regime -see below -depends on the temperature, room temperature lifetime data cannot be directly correlated with momentum spectra taken at the liquid nitrogen temperature.

New lifetime spectra have thus been collected simultaneously with momentum spectra by using a plastic (Pilot U, start) and a BaF 2 (stop) scintillator aligned at 90° from the axis of the Ge detectors, with the sample mounted on the cold finger of a cryostat at the centre of the cross. The resolution of the lifetime spectrometer was 257 ps, and the count rate in this case was 5 s -1 . The number of counts was 4×10 5 per spectrum. All spectra were analysed using the POSITRONFIT program [START_REF] Kirkegaard | PATFIT-88[END_REF] in one component after source subtraction. Since the spectra are expected to contain two or more unresolved components, the result of one-component fits must be intended as an effective average lifetime.

Coincidence Doppler broadening (CDB) of the annihilation radiation

All measurements were taken at the liquid nitrogen temperature simultaneously with lifetime spectra, as explained above, by means of two hyperpure Ge gamma detectors coupled to a multiparametric pulse analyser. Background rejection was achieved by requiring time coincidence within a window of 300 ns and fulfilment of the energy The momentum resolution was 3.7×10 -3 m 0 c. Twenty millions counts were collected in each spectrum. The analysis of the momentum spectra ( ) L p ρ was carried out, as explained in Refs. [16,17], by fitting a linear combination of momentum spectra measured in the same experimental condition for pure elements, as given by the following function:

conservation condition |E 1 + E 2 -2m 0 c 2 | < 2.
( ) V Zr Zr V Nd Nd V Y Y V Mg Mg bulk Mg ρ ρ ρ ρ ρ ρ C C C C F ) F 1 ( + + + + - = , (1)
where to the spectrum of Figure 2. Note that the vertical scale of the upper frame in Figure 3 is amplified by an order of magnitude for helping the visibility of the small signal given by trapping at vacancies in Mg. In order to eliminate uncertainties in the analysis of CDB spectra, a two-step fitting procedure was used. In the first step, free fits were used to determine the value of the trapping fraction F and its error as an average over a series of different fits giving approximately the same χ 2 . In the second step, the F value was kept as a constraint.

Figures 2 and 3 about here

Experimental results and discussion

Positron annihilation lifetime results

The evolution of the positron lifetime during ageing is shown in Figure 4 for the cold worked samples. For the purpose of comparison, results obtained from undeformed samples are also included in this figure. The horizontal axis of Figure 4 represents the ageing time normalised at the peak hardness ageing time (see arrows in Figure 1) which is equal to 22 h and 42 min, 5 h and 37 min and 3 h and 20 min respectively for the 0%, 6% and 12% deformations). The scaling of the time axis of Figure 4 corroborates that the deformation treatment affects the kinetics of the decomposition process without noticeable changes in the sequence of different stages of the precipitation process. Full symbols represent lifetime data taken at the liquid nitrogen temperature simultaneously with the CDB measurements. Open symbols are the lifetime data taken at room temperature for the same samples used by Macchi et al. low temperature [20] enables a larger fraction of positrons to reach the traps. The lowering of the positron lifetime at long ageing times, which is more evident at room temperature than at the liquid nitrogen temperature, certifies that the average distance between positron trapping sites is of the order of the positron diffusion length, i.e.

typically 100 nm at room temperature, and tends to increase with ageing time.

Figure 4 about here

Coincidence Doppler Broadening results

The momentum spectra ( ) p L = 0. The statistical noise has been reduced by averaging groups of three subsequent points from 1 to 2.6 atomic momentum units and of six points above 2.6 atomic units.

The solid lines through the experimental data are the best-fit curves obtained in accordance with the procedure described in Sec. 2.3.

Figure 6 about here

In the undeformed samples, the positron trapping is very weak (3%) after quenching and remains so even after 15 min at 250°C. On the contrary, the cold work increases the initial fraction of trapped positron by about 60% in the less deformed (6%) samples and near 70% in the more severely deformed (12%) samples. However, a very short heating leads to a sharp drop of the positron trapping. For example, for the 12% deformed samples the trapping fraction falls to about 15% after only 90 s at 250°C. After this initial decrease in positron trapping, further evolution of precipitation towards peak ageing entails an apparent increase in F by about 90%, which is then followed by a slight decrease in the case of the cold worked samples. As discussed above for the lifetimes, the decrease in positron trapping at long ageing time is due to coarsening of the precipitate distribution. As this effect is not seen in the undeformed sample, one can conclude that in this case, where the evolution is slower than that in the case of cold work, the precipitates are still dense at the scale of the positron diffusion length after 4 days at 250°C. However, this indication is not consistent with the lifetime data; thus the small rise of F at long ageing time might come from an imperfect fitting.

Additional information on the positron traps is given by the chemical composition of their environment. Figure 7 shows the total solute concentration probed by the trapped positrons, as given by the sum of the fitting parameters

Zr Nd Y sol C C C C + + =
, as a function of the ageing time scaled in proportion to the hardness peak time. This figure demonstrates that the intense trapping observed at the beginning of the ageing treatment is dominantly due to vacancy-like defects embedded in the Mg-rich matrix, without relevant decoration by solute ( sol C ≈ 0); this is opposite to the case of undeformed samples, where trapping is very weak and is due to vacancies deeply embedded in solute ( sol C near to 100%).

Figure 7 about here

Heating for 15 min does not substantially change both the local solute concentration and the trapping fraction in the undeformed sample. On the contrary, in the 12% deformed sample a short ageing of 90 s at 250°C rises the local solute concentration up to above 50%, in concomitance with a drop in the trapping fraction at its minimum of 15%. A jump is observed also for the less deformed (6%) sample after 15 min ageing (this sample was not tested at 90 s). It may be argued that the heating leads to rapid migration and loss of a large fraction of the vacancies that are introduced by the cold work. There is, however, a non-negligible fraction that survives and appears to be bound to solute clusters in an environment that is even richer in solute than expected for the stoichiometric compositions of the β', β Continuing the heat treatment makes the solute fraction gradually converge in all cases to an asymptote just below 20%. The asymptote is reached more or less in coincidence with peak hardening, when β', β 1 and β precipitates are densely distributed in the material as shown in the TEM images displayed in Figure 8.

Figure 8 about here

The information regarding the concentration of each solute species at the annihilation site is the following. Zr is positively identified only in the undeformed sample after quenching or after short ageing (15 min), when trapping is very small (F= 3%). In accordance with the interpretation given in Ref. [17], this result indicates that, in the absence of competition by other traps, a small fraction of positrons can reach and become trapped in sparse Zr dispersoids. The interesting information that can be drawn from this finding is that, in the undeformed material, after quenching and some dwelling time at room temperature, there are almost no free vacancies or vacancy-solute clusters containing Y and Nd, which in the model adopted in the present analysis, represents approximately all rare earths. On the contrary, these elements become well recognisable when the intensification of trapping indicates the beginning of precipitation. Figure 9 depicts the Y fraction

Zr Nd Y Y C C C C + +
. These data indicate that, excluding the pre-ageing conditions, the relative proportions of Y and Nd stay approximately constant in a ratio between 1:1 and 3:2 during the entire development of the ageing process (the horizontal line in Figure 9 is drawn at 58% that is the average value of the Y fraction for all samples, excluding the initial points).

Yttrium fracrtion (%)

Ageing tim e/Hardness peak tim e

  interactions during ageing of a cold worked Mg-RE based alloy F. Moia a , R. Ferragut a § , A. Calloni a , A. Dupasquier a , C.E. Macchi b , A.Somoza c and J. F. Nie d

  [14]). Coincidence Doppler broadening (CDB) of the positron annihilation radiation is another positron-based technique which can sensitively detect the identity of open volume defects and it is specially suited for identifying the nature of vacancysolute aggregates. The CDB technique is based on the measurement of the shape of the positron annihilation line at 511 keV with a very high signal/noise ratio (≈10 5 ), which leads to the determination of the momentum distribution of the electrons taking part in the positron annihilation. The shape of the electronic momentum contains

  1 keV, where E 1 and E 2 are the energies two detectors respectively. The component along the axis joining the two detectors of the momentum of the annihilating pairs is given by c

  the fitting parameters are the trapping fraction F (relative number of positrons that are trapped and annihilated in open volume defects) and the fractional concentrations C Mg , C Y , C Nd , and C Zr .Here it should be noted that the above fractional concentrations are to be interpreted as representing (with a possible distortion, due to the different positron affinity for the different elements[START_REF] Calloni | [END_REF] and to the lattice relaxation) the atomic composition in regions that are in immediate contact with the vacancy where the positron is trapped. It must be emphasised that the local composition in contact with a defect can differ from the average composition of the vacancy-solute aggregates. In Eq. 1, bulk Mg ρ is the experimental momentum distribution for bulk Mg measured from an annealed Mg with 99.9 % purity; distributions expected for annihilation in vacancies in pure elements, as obtained from measurements on cold worked samples after subtracting the contribution of annihilations from free positrons[16]. In our approximated analysis the contribution of rare-earths mixed to Nd is assumed to be included in the Nd term. The shape of the Mg bulk spectrum bulk Mg ρ is shown in Figure2. The

ρ

  for saturated positron trapping at vacancies ( X = Mg, Y, Zr, and Nd) are reported in Figure3as relative differences

  . The two sets of data are almost coinciding in the as-quenched condition and at the beginning of the heat treatment, but at later ageing stages the low temperature data are systematically above the room temperature data. This difference is consistent with positron trapping in a coarse distribution of extended defects like incoherent precipitates, which usually have open spaces at their interfaces with the magnesium matrix. In these conditions, referenced in the positron annihilation literature as diffusion-limited regime, the increase in the positron diffusion constant occurring at

  spectrum measured for bulk Mg shown in Figure2. This representation is used to enhance the details of the spectra in the high-momentum region, which is most important for the identification of the chemical species in contact with vacancy-like defects. Note that, for a fixed chemical composition of the environment of the positron traps, the curves scale vertically in proportion to the trapping fraction F. The symmetrical experimental curves have been folded around

Figure 5

 5 Figure 5 about here

Figure 9 Figure 1 :

 91 Figure 9 about here

Figure 2 .

 2 Figure 2. Momentum distribution spectrum for annealed Mg.

Figure 3 .

 3 Figure 3. Momentum distribution spectra for pure Mg, Y, Nd, Zr in conditions of high positron trapping at vacancies (relative difference to the spectrum of Figure 2 corresponding to annealed Mg).

Figure 4 .

 4 Figure 4. Mean positron lifetime as a function of the normalised ageing time (i.e., scaled in proportion to the hardness peak time) in cold worked and undeformed samples. Full symbols: room temperature (RT) measurements (after Ref. [11]). Open symbols: liquid nitrogen temperature (LNT) measurements (this work). The solid lines are for visual guide.

Figure 5 .

 5 Figure 5. Annihilation radiation spectra (relative difference to bulk Mg) in cold worked and undeformed samples at various ageing times. Solid lines: linear combination model.

Figure 6 .

 6 Figure 6. Comparison between the average positron lifetime (full symbols) and the trapping fraction F (open symbols) deduced from CDB data obtained from cold worked and undeformed samples vs. the ageing time. All measurements were performed at liquid nitrogen temperature.

Figure 7 .

 7 Figure 7. Solute fraction probed by the trapped positrons at the annihilation site in deformed and cold worked samples vs. the normalised ageing time. The solid lines are for visual guide.

Figure 8 .

 8 Figure 8. Transmission electron micrographs showing microstructuress typical of (a) 6% deformed samples aged for 8 hours at 250°C and (b) 12% deformed samples aged at 250°C for 4 hours.

Figure 9 .

 9 Figure 9. Ratio between Y and total solute atomic concentrations probed by the trapped positrons in cold worked and undeformed samples vs. the normalised ageing time.

Figure 1 F

 1 Figure 1

Figure 2 Page

 2 Figure 2

Figure 4 Page

 4 Figure 4

Figure 9

 9 Figure 9

  phases. It may be concluded that only the vacancies that are not stabilised by clustering with solute atoms are rapidly removed by heating.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	1 and β

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters