Stefan Sandfeld
email: stefan.sandfeld@kit.edu

Thomas Hochrainer

Peter Gumbsch

Michael Zaiser

Numerical Implementation of a 3D Continuum Theory of Dislocation Dynamics and Application to Microbending

Keywords: constitutive equations, dislocation dynamics, dislocation theory, plasticity of crystals user supplied): constitutive equations, dislocation dynamics, dislocation theory dislocation dynamics, continuum theory of dislocation, dislocation density, microbending

Crystal plasticity is governed by the motion of lattice dislocations. Although continuum theories of static dislocation assemblies date back to the 1950's, the line-like character of these defects posed serious problems for the development of a continuum theory of plasticity which is based on the averaged dynamics of dislocation systems. Only recently the geometrical problem of performing meaningful averages over systems of moving, oriented lines has been solved. Such averaging leads to the definition of a dislocation density tensor of second order along with its evolution equation. This tensor can be envisaged as the analogue of the classical dislocation density tensor in an extended space which includes the line orientation as an independent variable. In the current work we discuss the numerical implementation of a continuum theory of dislocation evolution that is based on this dislocation density measure and apply this to some simple benchmark problems as well as to plane-strain microbending.

Introduction

In recent years, the growing demand for physically motivated continuum descriptions of plasticity has led to a renewed effort to formulate continuum theories of dislocation kinematics and dynamics. Already in the 1950s, Kröner [START_REF] Kröner | Kontinuumstheorie der Versetzungen und Eigenspannungen[END_REF], Nye [START_REF] Nye | [END_REF], Bilby and co-workers [3] and Kondo [START_REF] Kondo | Proc. 2. Japan Nat. Congress of Appl. Mech[END_REF] independently formulated the classical continuum theory of dislocations based upon the definition of a second-rank dislocation density tensor (in the following denoted as Kröner-Nye tensor). The authors were well aware that the Kröner-Nye tensor, if used as a measure of the average dislocation state of a crystal, leads to a description of plastic deformation processes that is intrinsically incomplete. The main problem with classical dislocation density measures such as the Kröner-Nye tensor is that its averaging over a volume element containing many dislocation lines produces a measure that is invariant upon shear deformations of the averaging volume. As a consequence, any spatially homogeneous shear deformation can in principle not be reflected by the dynamic evolution of the Kröner-Nye tensor or similar measures. This renders the classical dislocation density measures highly problematic as foundations for a continuum theory of plasticity.

Nevertheless, plasticity theories have been formulated which are explicitly or implicitly based upon the classical dislocation density measures. As we shall discuss in detail, such theories, as for instance formulated by Acharya and co-workers [START_REF] Acharya | [END_REF], [6], [7] and by Sedláček et al. [8], apply to situations where dislocations are forming smooth bundles, i.e. adjacent dislocation lines are approximately parallel and have the same orientation. In more general situations, using the classical dislocation density tensor as the basis of a plasticity theory leads to non-closed formulations which need to be 'patched up' by phenomenological assumptions. Models which focus on the evolution of scalar densities of so-called statistically stored and geometrically neccessary dislocations (see e.g. [9] and references therein) are also essentially phenomenological. Such models can be considered as heuristic approaches towards formulating constitutive equations which allude to the physical mechanisms of dislocation microstructure evolution. However, the constitutive equations are not related to the evolution of a discrete dislocation system through systematic and mathematically meaningful averaging procedures, and the models therefore remain on the level of phenomenological expedients [10].

From the late 1990s on, Groma and co-workers developed a statistical approach towards deriving not only the kinematic evolution of dislocation systems, but also their internal interactions and stress driven dynamics, from systematic averages over ensembles of discrete dislocation systems [11], [12], [13]. This approach draws heavily on averaging techniques used in the statistical mechanics of interacting particle systems. It is therefore restricted to simplified systems of straight edge dislocations, which can be imagined as two-dimensional systems of point particles in the intersecting plane. In an average representation, these systems can be described in terms of densities of charged particles ('positive' and 'negative' dislocations distinguished by the two possible orientations of the line direction with respect to the Burgers vector). It is not obvious how to generalise these methods to threedimensional systems of connected dislocation lines, since the line-like nature of the dislocations and the continuous distribution of dislocation line directions need to be accounted for.

A description of dislocation systems in terms of higher-dimensional density measures dates back to Kosevich [START_REF]Kosevich Dislocations in solids[END_REF], who defines dislocation densities in a space which includes parameters characterising the line orientation as independent variable. In the context of dislocation dynamics, a similar approach was used in El-Azab's pioneering work [15], who extended the methods of statistical mechanics to systems of curved dislocation lines by using densities that evolve in a higher dimensional state space. However, while taking into account the orientation of line segments, both Kosevich and El-Azab do not directly account for the local connectivity of the segments. To illustrate the problem, we observe that their measures are not able to distinguish between a distribution of circular loops and a distribution of randomly oriented straight dislocation lines in a plane. The kinematic properties of these two systems are, however, very different.

Only recently, the mathematical foundations required for transferring the methods of statistical mechanics consistently to three-dimensional systems of curved dislocations have been formulated by Hochrainer's 'Extended Continuum Theory' (ECT) [16]. This theory relies on a geometrical description of dislocation lines and their averages, i.e. line densities, which uses the methods and formalisms of differential geometry. The theory is a direct generalization of Kröner's classical continuum theory, which is contained as a special case.

In this paper, we explore ECT from the point of view of its numerical implementation for a single slip geometry. Section 2 presents the theoretical background and is divided into three parts: section 2.1 briefly explains the constitutive relations governing the stress-strain response of a (macroscopic) solid. For completeness, we then introduce in section 2.2 Kröner's classical continuum theory and afterwards provide an overview of some of the theoretical aspects and basic notations of ECT in section 2.3. After discussing basic ideas of the numerical implementation in section 3, we explore in section 4 the evolution equations numerically and apply them to a continuous distribution of dislocation loops and, as a benchmark for kinematic consistency, to the expansion of a quasi-isolated single dislocation loop.

In section 5 we then apply the theory to analyze the dislocation evolution during bending of a free-standing thin film and investigate the size effects associated with a simple Taylor-type hardening model. Conclusions are presented in section 6.

Throughout the paper we use the following definitions and conventions: vectors and tensors are denoted by bold-face letters, whereas scalar quantities are written non-bold (e.g. α as opposed to α). Partial derivatives are abbreviated as ∂ x (•) := ∂(•) ∂x , ∇ denotes the gradient operator and the divergence of a vector (tensor) field is written as div; curl denotes the curl-operator. The vector product is denoted by ×. For double indices we assume the summation convention, if not stated otherwise. In c = a ⊗ b the tensor product is denoted by ⊗; this operation reads for cartesian coordinates in index notation c ij = a i b j . The symmetric part of a second rank tensor a can be obtained by sym a := 1 2 a + a T with a T denoting the transpose of a matrix.

Theoretical Background

Small strain elasto-plasticity

In continuum mechanics the deformation of a body is usually characterised by a smooth vector field of displacements u. The gradient ∇u is called the distortion tensor and is in the following denoted by β:

β := ∇u. (1)
The symmetric part of the distortion tensor is called the strain tensor:

ε := sym β = sym ∇u. (2)
In small-strain elasto-plasticity the distortion tensor is additively decomposed into a stress-free or plastic part β pl and an elastic part β el . Accordingly, the strain is additively decomposed as

ε = ε el + ε pl = sym β el + sym β pl . (3)
The local stress σ is related to the elastic part of the strain tensor by the consitutive law

σ = C : ε el = C : (ε -ε pl ), (4)
where C denotes the fourth order tensor of elastic moduli and ':' the doubly contracted tensor product, in Cartesian coordinates σ ij = C ijkl ε el kl . The last equation is accompanied by balance equations for momentum and angular momentum, which in the absence of body and inertial forces read

div σ = 0 and σ = σ T , (5) respectively.
The purpose of constitutive modelling in plasticity theory is to derive evolution equations for the plastic distortion β pl based on the current stress state and possibly various internal variables of the material. ECT is an example of a theory where the evolution of the distortion β pl is expressed in terms of internal variables which represent statistical averages over the discrete dislocation pattern within the crystal.

Classical continuum theory of dislocations

Kröner [START_REF] Kröner | Kontinuumstheorie der Versetzungen und Eigenspannungen[END_REF] defined the dislocation density tensor α as

α = curl β pl . (6)
A different but equivalent definition was given by Nye [START_REF] Nye | [END_REF]. In Eq. ( 6), the curl operator acts on the first indices of a tensor in the same way as it acts on a vector. From the definition of the Kröner-Nye tensor it immediately follows that

div α = 0, (7)
which reflects the physical fact that dislocations cannot end inside a crystal. When analyzing the information about the dislocation system that is contained in the Kröner-Nye tensor, it is important to specify the spatial resolution with which the plastic distortion is described. This is not always recognized in the literature, as the very concept of 'spatial resolution' seems odd when we are talking about a continuum theory. However, in order to understand the strengths and limitations of the classical continuum theory of dislocations it is indispensable to analyze how the information content expressed by the fundamental dislocation density measure changes as we move across scales. We start this analysis on the smallest scale on which a continuum theory can be used, i.e., we use continuous fields to describe the distortion of the crystal but assume that these fields are known with sufficient resolution such that individual dislocations can be 'seen' as singular lines.

Dislocation density tensor for a discrete dislocation system

If the spatial resolution of the theory is such that all dislocation lines are captured separately by the curl operation, the Kröner-Nye tensor completely characterises the dislocation system. In this case α can be explicitly related to the configuration of the discrete dislocation lines. We envisage a dislocation line as an oriented curve c(s) which we assume to be parametrised by arc length s, i.e. dc/ds is the unit tangent vector to the dislocation line. Furthermore we assume that all dislocations share the same Burgers vector b. In situations where dislocations of several slip systems are present, the following considerations hold for each slip system separately, and the total dislocation density tensor is obtained by summing over all slip systems. We formally define a density measure along the dislocation line by

δ c (r) = Lc 0 δ(c(s) -r)ds, (8)
where L c is the total length of curve c and δ(r) is the standard Dirac measure ('delta function') in three-dimensional space. Using this measure, we may write the discrete (indicated by superscript 'd') Kröner-Nye tensor as

α d = c δ c dc ds ⊗ b, (9)
where the sum runs over all dislocation lines in the system. The local rate of plastic distortion is given by Orowan's relation which in tensorial form reads

∂ t β pl,d = - c δ c v n ⊗ b = -v × α d . (10)
Here, v = vν is the local dislocation velocity, with v being the scalar velocity, and ν the unit vector in the dislocation glide direction. In writing down (10), we used that the glide plane normal is given by n = ν × dc ds . Defining the discrete dislocation current by J d = v × α d and using the definition of the Kröner-Nye tensor yields a kinematically closed evolution equation 1 for α d :

∂ t α d = -curl J d = -curl v × α d . (11)
We now attempt to use these relations for constructing a coarse-grained theory.

In doing so, we ask whether or not (and under which circumstances) direct coarse graining of the classical dislocation density measure leads to a loss of essential information.

Statistical averaging

We denote by (. . . ) V,r := (1/V ) Vr (. . . )d 3 r the spatial average over some small volume V r of size V centered at r (for convenience of notation, the subscripts V and r will be dropped in the following). A scalar measure of the average dislocation density is given by the dislocation line length within V r , divided by the averaging volume:

ρ = 1 V c c ∩V 1ds = 1 V c V Lc 0 δ(c(s) -r )d 3 r ds = c δ c . ( 12
)
1 We speak of a kinematically closed equation as we assume the dislocation velocity v at this point to be a given quantity. In general, v is a function of the local stress, which in turn depends on the dislocation arrangement. Hence, a mathematically closed theory requires additional relationships between the dislocation state, as expressed by α, and the dislocation velocity v. If v is a function of stress and line direction only, these relationships may be derived from Kröner's theory of eigenstresses [START_REF] Kröner | Kontinuumstheorie der Versetzungen und Eigenspannungen[END_REF]. This measure characterises the total dislocation density. Similarly, we define the average dislocation density tensor by

F
α = α d = c δ c dc ds ⊗ b , (13)
The average line direction is given by the unit vector

l = c δ c dc ds c δ c dc ds , (14)
and we define the geometrically necessary dislocation density ρ G by

ρ G = c δ c dc ds . ( 15
)
The ratio ρ G /ρ is always smaller equal than 1 because of the triangular inequality.

The average dislocation density tensor can be written in terms of ρ G and l as

α = ρ G l ⊗ b. (16)
In general, the average dislocation density tensor (16) does not fulfill the equation

∂ t β pl = -v × α, (17)
if we interpret v as the pointwise vectorial average of the velocities along the dislocation lines. By use of (10) together with averaging of (11), we can attempt to obtain an evolution equation for the averaged tensor. Writing this out, we obtain

∂ t α = -curl c δ c v × dc ds ⊗ b . ( 18
)
The threefold product within the average can be written as a product of averages if, and only if, all dislocations within the averaging volume share the same tangent vector l = dc/ds and velocity v. This is possible either if only one dislocation is present (the discrete case) or if the dislocations form smooth line bundles. In these cases, ρ = ρ G , α = ρ l ⊗ b, and (11) holds both on the local and on the averaged scale.

In the general case, however, the averaging volume contains dislocations of different orientations. Thus, averaging leads to a reduced dislocation density ρ G < ρ, and the dislocation density tensor does not obey (11) with v understood as the average velocity. Instead, additional terms appear in the evolution equation. The necessity of accounting for such terms, which have the formal structure of correlators between the dislocation density and velocity/line direction, has been recognized in the literature (see e.g. [17]). However, the derivation and closure of a higher-order theory which properly accounts for such terms requires a formidable theoretical effort. Until now, published attempts either remain at the level of declarations-of-intention, or provide phenomenological patches which do not really resolve the underlying theoretical problems. To resolve the problem discussed above, the ECT distinguishes dislocation lines in a given spatial point r according to their line direction l. As long as we are on the level of a discrete description, this is simply redundant. However, the physical rationale for this distinction becomes evident when we average over a mesoscopic volume: While it is in general unrealistic to assume that all dislocations contained in a mesoscopic volume have the same direction and move with the same velocity in response to an acting stress, it is much more realistic to assume that those dislocations which do have the same direction move in the same manner.

F

If dislocations move by glide only, their line direction can be parametrised by a scalar variable, for instance the angle ϕ between line direction and Burgers vector. A point in configuration space is thus denoted by (r, ϕ). With regard to the numerical examples treated subsequently, we restrict ourselves to introducing ECT only for this situation. We note, however, that the theory is not restricted to this case but can be generalized to include out-of-plane dislocation motion and multiple slip systems (cf. [16]).

Definition of the second order dislocation density tensor

To define the dislocation density tensor and its averages in the three-dimensional configuration space we proceed in close analogy to the classical theory. Unlike in the classical case, however, we consider averages over so-called lifts of dislocation lines to the configuration space instead of averages over the spatial dislocation lines themselves. In order to introduce the concept of lifted curves we define a coordinate system such that the 1-direction points into the direction of the Burgers vector and the glide plane is spanned by the 1-and 2-directions. To each point c(s) = (c 1 (s), c 2 (s)) of a dislocation line we assign the angle ϕ between the tangent dc/ds and the Burgers vector, that is ϕ(s) = arctan( dc 2 ds / dc 1 ds ). We define the lift C of a given curve c to the three dimensional configuration space as

C(s) = (C 1 (s), C 2 (s), C 3 (s)) := (c 1 (s), c 2 (s), ϕ(s)). (19)
Figure 1 gives a visualisation of this concept. We note that the tangent to the lifted curve

dC ds = dc 1 ds , dc 2 ds , dϕ ds = dc ds , k(s) , (20)
contains as third coordinate the curvature k(s) of the curve at c(s). As a consequence of this implicit definition of a metric the lifted curves C are not parametrised by arc length in the configuration space. 1 By transferring the definition of δ c to the lifted curve C as

δ C (r) = C δ(C(s) -(r, ϕ))ds, ( 21
)
1 It is worth noting that the definition of the second order dislocation density tensor does not necessarily require a metric or volume element. However, an invariant definition requires the use of advanced mathematical concepts as e.g. differential forms. We refrain from introducing these concepts in the current paper and refer the interested reader to [16] and especially [START_REF] Hochrainer | Evolving systems of curved dislocations: Mathematical foundations of a statistical theory[END_REF] for a more thorough treatment. we define the dislocation density tensor of second order in complete analogy to the classical one as

α II (r,ϕ) = C δ C dC ds ⊗ b . ( 22
)
The sum is again taken over all dislocations in the system, and the averages are taken over a volume element in configuration space which is centered at r and ϕ.

It is illustrative to express also α II as the product of a density function ρ (r,ϕ) and a generalised tangent which we denote by L (r,ϕ) . The definition of these two objects is not completely analoguous to the classical ones. The density function is defined as

ρ (r,ϕ) = C δ C dC ds r , (23)
where • r measures the length of the spatial projection of a vector in the configuration space. The tangents to the lifted curves have the property

dC ds r = dc ds = 1. (24)
In ( 23), the averaging volume -a volume element in configuration space -contains dislocations of one orientation only. Therefore, there is no cancellation -e.g. of dislocations of opposite directions -during averaging of the spatial tangent vectors. Therefore, the norm • r may be interchanged with the averaging. We thus find that ρ (r,ϕ) gives the spatial line length per volume (of the configuration space) of dislocations at r with direction l(ϕ) = (cos ϕ, sin ϕ). The generalised tangent is defined as

L (r,ϕ) = C δ C dC ds ρ (r,ϕ) . ( 25
)
We note that the first two components of L (r,ϕ) are just the canonical spatial direction at ϕ and the third component contains the average curvature k (r,ϕ) , that is

L (r,ϕ) = (cos ϕ, sin ϕ, k (r,ϕ) ) = (l(ϕ), k (r,ϕ) ). ( 26
)
From the definitions ( 22), ( 23) and ( 25) we easily find

α II (r,ϕ) = ρ (r,ϕ) L (r,ϕ) ⊗ b. (27)
The classical dislcoation density tensor α can be evaluated from the density function ρ (r,ϕ) as The second order dislocation density tensor may be represented by a 3×2 matrix which reads

α (r,ϕ) = 2π 0 ρ (r,ϕ) l ϕ dϕ ⊗ b. (28)
α II (r,ϕ) = ρ (r,ϕ)    b 1 cos ϕ b 2 cos ϕ b 1 sin ϕ b 2 sin ϕ b 1 k (r,ϕ) b 2 k (r,ϕ)    = ρ (r,ϕ)    cos ϕ 0 sin ϕ 0 k (r,ϕ) 0    b, (29)
where the latter follows from the chosen coordinate system for which we took the 1-direction parallel to the Burgers vector and hence find b = (b, 0). We remark that in analogy to (7) also the generalised divergence of α II must vanish (cf. [16]), which likewise reflects the physical fact that dislocation lines do not start or end inside a crystal. With (29) this condition reads cos ϕ∂ x ρ + sin ϕ∂ y ρ + ∂ ϕ (ρk) = 0.

(30)

The generalised velocity

In the preceding section we defined the dislocation density tensor α II as an average of density functions characterising discrete lifted curves. In order to obtain the evolution of α II , we need to consider the velocity of the lifted curve in the configuration space. As may be seen in Figure 1 the velocity of the lifted curve contains, besides the spatial (horizontal) velocity v orthogonal to the spatial dislocation line, also a component in the ϕ direction which accounts for the rotation of line segments during dislocation motion.

It is straightforward to show that the pseudo-scalar rotation velocity ϑ of a moving curve parametrised by arc length s is determined by the change of velocity along the line. For the pseudo-scalars v, k and ϑ we use the following sign convention:

We consider the pseudo-scalar velocity of a positively oriented loop as positive if the loop expands (compare Fig. 2). The curvature of such a loop is also considered as positive. From this sign convention it follows that the pseudo-scalar rotational velocity is

ϑ(s) = - dv(s) ds . (31)
Using these conventions, we define the generalised velocity at a point C(s) of a lifted discrete curve as

V d (s) = v sin ϕ(s), -v cos ϕ(s), ϑ(s) . (32)
We note that unlike the spatial velocity, which is orthogonal to the dislocation line, the generalised velocity is in general not orthogonal to the lifted curve. The evolution of α II

In analogy with the classical definition we define the generalised (discrete) dislocation current

J II,d = C V d × δ C dC ds ⊗ b. ( 33
)
The evolution of the averaged dislocation density tensor of second order is obtained in analogy to [START_REF] Hochrainer | Evolving systems of curved dislocations: Mathematical foundations of a statistical theory[END_REF]:

∂ t α II (r,ϕ) = -curl J II,d (r,ϕ) . ( 34
)
The curl operator again acts on the columns of the coefficient matrix of a tensor. The difference to the evolution of the classical dislocation density tensor is that in (34) the product appearing in the generalised dislocation current (compare (33)) can be understood as a product of averages under much weaker assumptions on the dislocation configuration. The most important assumption in this respect is that nearby dislocations of the same line direction have the same curvature, which is reasonable in quasi-static situations, where the curvature essentially balances the local stress. Under this assumption, the average dislocation current can be written as where the average velocity V (r,ϕ) = V d can be expressed in terms of the average pseudo-scalar velocity v (r,ϕ) = v and rotation velocity ϑ (r,ϕ) = ϑ as

J II = V (r,ϕ) × α II (r,ϕ) , (35)
V (r,ϕ) = (v (r,ϕ) sin ϕ, -v (r,ϕ) cos ϕ, ϑ (r,ϕ) ) . (36)
Consequently, the evolution equation for the dislocation density tensor of second order is obtained as

∂ t α II (r,ϕ) = -curl V (r,ϕ) × α II (r,ϕ) . (37)

Evolution equations for the scalar field variables

As we saw in (29), in the single glide situation under consideration α II is completely defined by the two scalar fields ρ (r,ϕ) and k (r,ϕ) . Hence, the evolution equation for α II can be translated into evolution equations for these functions. As the derivation of these evolution equations involves either abstract methods or simple but lengthy calculations we only give the result and refer the reader to [16] for details:

∂ t ρ = -div (ρv) + ∂ ϕ (ρϑ) + ρvk, (38)
∂ t k = -vk 2 -∇ 2 L (v) -∇ V (k). ( 39
)
Here, ∇ 2 L and ∇ V denote the second derivative along the generalised line direction and the first derivative along the generalised velocity direction, respectively. To avoid cumbersome notation we dropped the index (r,ϕ) . In the following if no indices occur we always will imply that all quantities are evaluated at the point (r, ϕ) in the configuration space R × R × [0, 2π).

Evolution equation for the plastic distortion

The evolution equation for the plastic distortion β pl takes a similar form as (10) in the classical formulation:

∂ t β pl = 2π 0 ρ (r,ϕ) v (r,ϕ) dϕ n ⊗ b. ( 40
)
This quantity provides the link to the material response on the continuum level as introduced in section 2.1.

Numerical Implementation

Geometry

As in the preceding section we assume that all dislocations are moving on a single glide system with glide plane normal n = e 3 and Burgers vector b ≡ be 1 . The system is assumed homogeneous in the e 3 direction, i.e., we consider a statistically homogeneous assembly of active glide planes with constant average spacing. We define a discrete time step ∆t := t i+1t i with i ∈ N + . A forward Euler scheme is used for numerical time integration of the evolution equations. Hence, only values at step t i are needed to compute the new values at step t i+1 ; no additional history variables need to be stored. This is beneficial for the required computer memory, which is an issue as the multidimensional structure of the configuration space requires a large amount of information to be stored for multiple fields. Furthermore, we use a fixed mesh which demands a relatively fine space discretisation with corresponding large storage requirements. The forward Euler time integration scheme is an unstable scheme which is only first order accurate and suffers from dispersion. Our main reason for using this scheme is that its simple explicit structure allows for an easy identification of the separate effects stemming from the different terms in the evolution equations. Implicit time integration schemes, on the other hand, require additional linearization steps along with the iterative solution of the resulting difference equations. It is then not straightforward to distinguish effects stemming e.g. from the discretisation, the initial values, the approximative nature of the Newton linearisation or the partial differential equations themselves. For the explicit Euler scheme used together with a compatible choice of numerical derivatives, we could identify a 'working regime' within which this integration method was numerically reliable and could be used to 'test' more fundamental aspects of the numerical implementation1 . The 3-dimensional configuration space is discretised by a uniform mesh. It turned out that a sufficient resolution in the angular direction requires between 60 and 120 nodes. In the case of quasi-discrete lines we use a discrete Gaussian distribution to approximate the Dirac delta function. Since we use a fine spatial resolution for this representation (e.g. about 60 nodes to discretise the part of the Gauss function which is > 0.1% of the peak value) dispersive effects can be neglected. Derivatives w.r.t. the configuration space are approximated by finite differences. Derivatives which govern transport were approximated by an upwind scheme. This is a numerical scheme which uses information about the flow direction (e.g. the direction of propagation of a wave) to determine whether to use a forward or backward difference stencil, such that only information from upstream of the flow is used. Similar schemes are widely used to solve hyperbolic partial differential equations in computational fluid dynamics and reduce oscillations in direction of wave travel (e.g. [START_REF] Harten | [END_REF]).

The upwind method based on first order forward and backward differences is oscillation-free but at the cost of a very strong 'smearing-out' effect in direction of wave travel caused by a large amount of 'numerical viscosity'. We therefore use a scheme based upon second order accurate forward and backward differences, which is of much better quality than the first order scheme but slightly oscillatory [START_REF] Strikwerda | Finite difference schemes and partial differential equations[END_REF]. As a very simple remedy against the undershoot effect causing the oscillations, we cut off the (unphysical) negative density values and distribute this negative amount over all nodes with positive density. This procedure causes a steepening gradient on the downstream side of the wavefront and a slight increase of the wave's maximum -both of which are only minor effects in the cases we study.

All derivatives for nonconvective (e.g. diffusion) terms were approximated by central difference schemes of second order accuracy. To determine the time step size we use the Courant-Friedrichs-Lewy condition [START_REF] Morton | Numerical solution of convection diffusion problems[END_REF]

v ∆t ∆x ≤ 1, (41)
which serves as a (necessary) stability criterion for purely advective equations of the form

∂ t ρ + v∂ x ρ = 0 (42)
that are solved by explicit time integration schemes such as the forward Euler scheme. In our simulations we accordingly choose ∆t ≈ 0.1∆x/v where v is the largest occurring velocity.

Test Cases

In this section we investigate two examples, both of which are simple enough such that we still can easily predict the outcome while exploring the evolution equations of ECT. The first example treats a homogeneous distribution of circular loops. This case, which can be solved analytically, illustrates some basic properties of the theory. The second example treats a single quasi-discrete circular dislocation loop. This serves to illustrate and address issues related to the numerical discretisation of the system in a higher-dimensional configuration space.

Homogeneous distribution of expanding loops

For a homogeneous distribution of equidistant loops with the same initial radius r 0 in an unbounded glide plane, all derivatives of ρ, k and v w.r.t. the configuration space vanish. Then the system of evolution equations ( 38) and (39) reduces to

∂ t ρ = ρvk, (43)
∂ t k = -vk 2 . ( 44
)
(43) only consists of a production term, which yields the change of line length during expansion of the loops, while (44) governs the change of curvature k(t) = 1/r(t) due to the expansion. Note that here ρ and v are independent of ϕ. Thus, we can directly obtain the evolution of the plastic distortion from (40) as

∂ t β pl = 2πρbv. ( 45
)
We could have obtained the same evolution equation for the curvature by considering the change of radius of a loop with initial radius r 0 which expands with the velocity v:

k = 1 r(t) = 1 k 0 + vt -1 (46)
⇒ ∂ t k = -v 1 k 0 + vt -2 = -vk 2 . ( 47
)
The evolution equation for the scalar density can then be obtained by considering homogeneously distributed loops with mean center-to-center spacing d and initial radius r 0 . The density is the line length 2πr(t) per volume, i.e. in terms of k ρ = 2π k

1 d 3 . ( 48
)
We obtain the evolution equation as time derivative of (48):

∂ t ρ = 2π d 3 ∂ t 1 k = -ρ ∂ t k k = ρvk. (49)
Hence, the case of a homogeneous distribution of equidistant, expanding or shrinking loops with the same initial radius is correctly represented by the theory.

In our introduction we mentioned that an essential shortcoming of Kosevich's and El-Azab's continuum theories is manifested by the fact that the kinematic behavior of a distribution of straight dislocation lines is not distinguished from that of a distribution of loops with the same radii. The ECT, however, reproduces the correct increase of line length for the evolving system of distributed loops, while for the case of the straight line distribution the theory yields pure advection without change in line length.

Quasi-discrete expanding circular loop

We again assume constant velocity and curvature such that there is no rotational velocity and no change of curvature along the line. As all derivatives of v and k are zero the evolution equations now simplify to

∂ t ρ = -div (ρv) + ρvk, (50)
∂ t k = -vk 2 . ( 51
)
The first part of the evolution equation for ρ governs spatial transport, whereas the second part accounts for changes in line length due to expansion or shrinkage of curved dislocation segments (here forming circular loops). The change of curvature which goes along with the expansion or shrinkage of segments is determined by the evolution equation for k.

Numerical results: The problem of line fragmentation

As an approximation to a discrete dislocation loop we generate a continuous density function by replacing the δ-distribution in [START_REF] Sedláček | [END_REF] with a Gaussian distribution. The spatial projection of the initial distribution, that is the total dislocation density, is depicted in Fig. 3(a). Time integration of the evolution equations ( 50) and (51) leads to a spatial density distribution as shown in Fig. 3 While the total amount of density integrated over the plane is -up to a numerical error -predicted correctly (c), the loop heavily fragments. This also displays in the evolution of the numerical divergence depicted in (d).

the total dislocation line length, as expressed by the volume integral of the density function over the configuration space, is described correctly (Fig. 3(c)), except for a small deviation which is caused by discretisation errors. Nonetheless, the density distribution gets spatially fragmented as the loop expands. This fragmentation is directly related to a growing divergence of α II . As the time derivative of α II is a curl (compare (37)) the evolution equations themselves conserve the solenoidality of α II . However, finite difference schemes in general can not guarantee that the discrete curl -and div -operators fulfill the relation div curl = 0. This shortcoming manifests itself in a strong increase of the mean divergence (1/N nodes ) div α II as shown in Fig. 3(d). This increase mirrors the visible fragmentation of the dislocation density distribution during expansion of the loop. This fragmentation is unphysical. Moreover, in the case of inhomogeneous velocity fields the temporal evolution depends on derivatives of the density along the generalised line direction. The numerical evaluation of these derivatives becomes completely inaccurate if loop fragmentation becomes pronounced. We note that other authors have observed similar problems, see e.g. [START_REF] Sedláček | [END_REF] where the problem is addressed by utilising polar coordinates. This avoids fragmentation, but the approach can be used only for the special case of loops with the coordinate origin as their center, and not for the evolution of general distributions. The remainder of this section deals with generic strategies that avoid fragmentation within the framework of the chosen finite difference scheme.

Direct relaxation of the dislocation density tensor α II

A straightforward method to maintain solenoidality of the density distribution is to minimize the divergence of α II by means of a sequence of 'relaxation' steps following each time step. This process ensures iteratively that div α II as a measure of fragmentation is kept below a given tolerance. We minimise the total divergence of α II using the iteration formula

α II new ⇐ α II old + λgrad div (α II ), ( 52
)
where λ is a factor controlling the step size, which can be adjusted to achieve efficient relaxation while avoiding 'overshoots'. However, this method has clear drawbacks: The operators grad div (•) introduce second derivatives which have a diffusive effect and suppress fragmentation by homogenizing the dislocation density distribution along the dislocation line. However, the diffusive effect also acts in the perpendicular direction and causes a broadening of the line. Furthermore, the iteration scheme (52) conserves α II and not ρ, causing unphysical changes in dislocation density during the relaxation step. An improved relaxation scheme which conserves the total density and avoids diffusive 'flattening' of the dislocation density distribution is obtained by modifying (52) such that ρ is relaxed only along the line:

ρ new ⇐ ρ old + λ∇ 2 l ρ, (53)
where the relaxation is again carried out to minimize div α II . Irrespective of the method used, relaxation of the system is computationally expensive. Achieving a well-relaxed configuration requires multiple relaxation steps, each of which is computationally about as expensive as one time integration step for the evolution equations. To obtain a computationally more efficient method for preserving solenoidality we take a closer look at the mechanism that leads to the fragmentation and devise a remedy which restores the unfragmented state without introducing an additional iterative process.

Correction of line fragmentation by tangential diffusion

To better understand loop fragmentation from a geometrical point of view, we consider a curved dislocation line in the configuration space (for simplicity we envisage a circular loop but the argument extends to any curved line). To obtain a continuous and differentiable density distribution we need to approximate the Dirac δ-function representing the discrete line in terms of continuous and differentiable functions. Doing this on a discrete grid amounts to a discrete convolution with an approximation function, i.e., we map (a subset of) R to (a subset of) N. For poor resolution in ϕ-direction and a better resolution in x and y direction, Fig. 4 shows the density distribution obtained with Gaussian functions. It is important to realize that due to the mapping, each of the 'blobs' representing segments of the dislocation line is located on a particular ϕ-plane. As all points of the blob move into the same direction v, which is perpendicular to the spatial line direction l(ϕ), the blobs are drifting apart during loop expansion and the line becomes fragmented. We consider the effect of loop expansion at two successive time steps t i and t i+1 and its impact on the change of arc length of a line segment s i for the case of a continuous lifted loop and for a lifted loop with discretized orientations. Fig. 5 shows the geometrical relations between these two cases. To understand how the differences affect dislocation density evolution, we may then think of a single lifted line as as a bundle of parallel lines with partial Burgers vectors and characterize the bundle by a space-dependent density ρ.

Let us first consider the expansion of a continuous lifted loop. As the loop expands in the spatial plane a circle segment s i with cone angle ∆ϕ gets elongated to s i+1 (cf. Fig. 5(a)). Expressing this relation in terms of the curvature we get

k i = ∆ϕ s i and k i+1 = ∆ϕ s i+1 . (54)
The increment of spatial arc length ∆s from step i to step i + 1 then is defined through

∆s := s i+1 -s i = ∆ϕ 1 k i+1 - 1 k i . ( 55
)
The arc length increase stems from the fact that the line orientation along the circle segment varies between ϕ 0 -∆ϕ 2 and ϕ 0 + ∆ϕ 2 , cf. Fig. 5(a), leading to divergent trajectories and a separation of the endpoints of the segment. For a bundle of lines, the local density within the bundle remains unchanged while the overall density increases as the bundle occupies more space.

The numerically discretised loop, on the other hand, consists of 'straight' line segments (the orientation is constant within each segment). Hence, each point of the segment translates into the same direction and no increase in line length ∆s can take place (thick black lines in Fig. 5(b)). However, ECT is based on the description of continuous lifted lines and therefore its equations are constructed to generate the correct additional amount of density (which is governed by the production term ρvk of (38)). Unfortunately, after discretization the transport terms fail to locate this amount in the right place.

The consequences are illustrated in Fig. 6: In case of a bundle of continuous lifted lines, Fig. 6 (a), additional line length is introduced by stretching the curved segments while the density remains constant. In case of a bundle with discrete orientations, Fig. 6 leading to fragmentation, while the additional line length causes an increase in density,

∆ρ = s i+1 s i -1 ρ = ∆s s i ρ. ( 56
)
This suggests to remedy the situation by redistributing the additional density along the line.

The intuitively obvious method to do this is shown in Fig. 6(c): We might remove the added density ∆ρ and distribute it to the left and right of the original segments, such as to extend these segments by exactly the amount needed to restore connectivity. Unfortunately, numerical implementation of this idea is not straightforward since we are working with a fixed spatial grid where the notation of 'segment length' is not defined. Instead, we implement a diffusive spreading scheme: We continue each segment along its direction into a sequence of segments of similar orientation and length as shown in Fig. 6(d), and distribute the excess density on its two nearest neighbors in this sequence. Conversely, if these neighbors already carry density (and thus produce additional density ∆ρ during loop expansion) half of the newly generated density from both neighbors is distributed back onto the original segment.

Thus, at each time step the density at r is diminished by the newly created density ∆ρ(r) but increased by the contributions ∆ρ L (r) and ∆ρ R (r) from the neighbors at r ± s i l,

∆ρ L (r) = ∆s 2s i ρ(r -s i l), ( 57
) ∆ρ R (r) = ∆s 2s i ρ(r + s i l). ( 58
)
The total density change due to this diffusive re-distribution is then given by ∆ρ diff (r) = ∆s 2s i ρ(r -

s i l) + ρ(r + s i l) -2ρ(r) . ( 59
)
Taylor expansion around r of first two density terms on the right-hand side yields

ρ(r -s i l) ≈ ρ(r) -s i ∇ l ρ(r) + s i 2 2 ∇ 2 l ρ(r) (60) ρ(r + s i l) ≈ ρ(r) + s i ∇ l ρ(r) + s i 2 2 ∇ 2 l ρ(r), (61)
where ∇ l and ∇ 2 l are the first and second derivatives along the spatial direction of the line. Inserting the above two equations into (59) yields

∂ t ρ diff (r) = lim ∆t→0 ∆s ∆t 1 2 s i ∇ 2 l ρ(r) , (62)
which is the general formulation of the diffusive correction term. ters. We use the identity

∆s s i = - ∆k k i = -∆t ∂ t k k i for ∆t → 0 (63)
and observe that the segment length s i is related to the curvature and orientation space discretisation parameters by s i = ∆ϕ/k i . This allows us to rewrite (62) as

∂ t ρ diff = -µ 2 ∆ϕ k 2 ∂ t k k ∇ 2 l ρ. ( 64
)
In this expression we have introduced µ as a shape factor to account for the fact that, while our argument is based on considering the lifted line as a sequence of seqments of a parallel line bundle, a numerical discretization of the line in terms of Gaussian functions is more akin to a sequence of rounded 'blobs'. In practice, we determine this factor from numerical experiments to achieve an optimum balance between the dual requirements that (a) fragmentation should be efficiently suppressed and (b) the density profile of the line bundle should be preserved, i.e., no additional dispersion in the glide direction should be introduced. For the chosen discretisation of 60 nodes in angular direction an optimal value for µ can be identified as

µ opt ≈ 0.3. ( 65
)
For larger values µ > µ opt the diffusive effect introduces a significant reduction in peak density and smearing out of the line profile, whereas for µ < µ opt the violation of div α II = 0 is still appreciable. In Fig. 7(b) one can observe the mending effect of the correction term. In Fig. 8 the temporal evolution of the sum of the mean absolute value of div α II is shown for different correction methods. Relaxation of the system along the dislocation line proves the most effective method, however, carrying out the relaxation to achieve the shown accuracy requires between 10 and 30 relaxation steps after each time step, which causes an increase of the overall computation time by a factor of more than 10. Evaluating the diffusive correction, on the other hand, can be done 'on the fly' and does not cause a notable increase in computation time. For a geometry factor of µ = 0.3, we observe an efficient supression of loop fragmentation.

We finally note that the same considerations as above also hold for the shrinkage of a loop. In this case the diffusion coefficient becomes negative which allows the loop to localise.

Application to Microbending

To investigate the performance of our numerical scheme in a physically relevant situation, we consider the problem of micro-bending of a free-standing thin film. This problem has several interesting features: (i) One is dealing with an intrinsically heterogeneous deformation state where the 'composition' of the dislocation arrangement changes over time due to the growing strain gradients; (ii) These strain gradients are associated with pronounced size effects that have been studied extensively using various dislocation-based models (see e.g. [9], [10], [START_REF] Sedláček | [END_REF], [26]); (iii) the presence of an elastic core at the centre of the specimen leads to the spontaneous emergence of large curvatures in the dislocation system, and hence provides Ideally, the sum would be zero for all steps. Due to the numerical approximation, however, it is non-zero right from the beginning and always increases. A larger amount of diffusion decreases the sum of div α II by comparison. The optimum amount of diffusion obtained for µopt is a trade-off between minimizing div α II and maintaining the correct peak density value. Relaxation of ρ yields the best result but at a very high computational cost.

an interesting testing case for the numerical accuracy of our computations in a situation where all terms in the evolution equations ( 38) and (39) yield appreciable contributions to the evolution of the dislocation system. We consider bending of a free standing thin film of thickness h. The normal vector of the free surfaces is denoted by n s , and we assume without loss of generality that n s is perpendicular to the y axis of a Cartesian coordinate system. The dimensions of the film in the directions perpendicular to n s are assumed large so that the system can be considered homogeneous in these directions. The film is deformed by bending around an axis parallel to the y direction. Plastic deformation by dislocation glide may occur on two slip systems that are symmetrically inclined with respect to n s as shown in Figure 9. The Burgers vector of each slip system is perpendicular to the y axis, defining a plane-strain geometry. Since both slip systems are equivalent, it is sufficient to analyse the evolution of the dislocation system for only one of them. Furthermore, we may restrict ourselves to the evolution of strain and dislocation density on a single slip plane since the system is homogeneous in the directions parallel to the film. We define a Cartesian coordinate system such that the slip plane coincides with the plane z = 0 as shown in Figure 10 and the origin is located on the neutral fiber in the center of the film. The Burgers vector of the considered slip system points in the positive x direction. The shear stress in the considered slip system can be formally envisaged as the sum of a reference shear stress τ 0 which describes the stress state in a film that is bent to the same radius of curvature R in a purely elastic mannner, and a shear stress τ 1 related to the plastic strain γ, which is the sum of the plastic strains from both slip systems. The stress τ 0 is given by

τ 0 (x) = sin(2θ)h R(1 -ν) G x l , ( 66
)
where G is the shear modulus of the material, ν is Poisson's number, and l is the film width projected on the x direction (Figure 10). For the geometry under consideration, plastic deformation of a fiber leads to a stress reduction that is proportional to the axial strain. Assuming isotropic material properties and making the standard assumption that straight specimen cross sections (cross sections parallel to the bending axis) remain straight during bending, we find that

τ 1 (x) = - sin 2 (2θ) 1 -ν Gγ(x), (67)
where we have accounted for the plastic strain contributions from both slip systems.

The bending moment M (moment per unit length in the y direction) corresponding to a given bending radius R and plastic deformation state is evaluated from the stresses τ 0 and τ 1 as

M h 2 = 4 sin(2θ) 1/2 0 (τ 0 (u) + τ 1 (u))u du, ( 68
)
where u = x/l and τ 0,1 (u) := τ 0,1 (x = ul). The scaled bending moment M/h 2 in the following serves as measure of the characteristic stress in our bending simulations.

Constitutive relations

We use a linear-viscous model of overdamped dislocation motion by assuming that the dislocation velocity depends linearly on the difference between the local shear stress τ and a local yield stress τ y :

v =      b B (τ -τ y ) if τ > +τ y , b B (τ + τ y ) if τ < -τ y , 0 otherwise.
(69)

Physically, the yield stress on a given slip system is governed by the intersection of moving dislocations with 'forest' dislocations of other systems, leading to junction formation. In the present case, both slip systems are equivalent, hence we can relate the yield stress to the dislocation density on the slip system itself. We assume a Taylor relationship,

τ y (x) = aGb ρ(x), (70)
where a ≈ 0.4 is a nondimensional constant characterising the strength of dislocation-forest interactions.

In our constitutive model we neglect effects of large-scale dislocation curvature on the dislocation velocity (line tension as considered in [START_REF] Sedláček | [END_REF]26]). We also neglect back stresses resulting from the 'piling up' of dislocations of the same sign [26,27], which are also implemented in DDD simulations of microbending ( [24], [25]). . Consistent with our assumption that dislocation motion is controlled by forest interactions, we do not impose specific boundary conditions to account for image effects at the free surfaces. Instead, we allow dislocations of all orientations to enter or leave freely, simply extrapolating the dynamics inside the film across the boundary (similar to the 'natural' outflow boundary conditions in [24], wherecontrary to our boundary conditions -only outflow and no inflow of dislocations from the surface was allowed). In physical terms our boundary conditions imply that we assume that surface sources have no activation stress other than the yield stress that controls near-surface dislocation motion inside the sample.

Results

As initial condition for our simulated bending tests, we assume an isotropic dislocation pattern with zero mean curvature and space-independent total dislocation density ρ = 10 13 m -2 , i.e., initially dislocations of all orientations are present with equal density ρ(x, ϕ) = (1/2π) × 10 13 m -2 . This initial condition describes a homogeneous and statistically isotropic arrangement of straight dislocation lines threading the film at random points. We carry out quasistatic simulated bending tests by increasing the load τ 0 (x) in small steps. In physical terms, this corresponds to slowly decreasing an imposed curvature radius on the thin film. We start from the critical stress level where the shear stress at the surface reaches the yield stress corresponding to the initial dislocation density. After each stress increment, we solve the equations of motion while simultaneously evaluating the increase of the plastic strain and the concomitant decrease of local stresses and increase in flow stress. Due to these changes in the internal stress state the strain rate gradually decreases towards zero. We trace this relaxation until the maximum strain rate has everywhere dropped below a prescribed low level and then record the bending moment M and strain profile γ(x) as well as the dislocation and curvature patterns The following examples were computed using finite differences with an equidistant mesh. The resolution was about 700 spatial nodes and 100 angular nodes; a fifth-order Runge-Kutta method with adaptive step size control turned out to be an effective time integration scheme: the computational cost on a 2.4 GHz quad core work station for each of the following simulations is about 20 minutes.

Typical density and curvature patterns are shown in Fig. 11. These patterns are characterised by the presence of an elastic core region where the stress τ 0 (x) is less than the yield stress corresponding to the initial dislocation density. No deformation activity takes place within this region and, hence, the dislocation density and curvature remain at their initial values. Curvature spatially localizes near the boundaries of the elastic core region which are moving inwards with increasing stress. In terms of the orientation coordinate, curvature is strongest for near-screw orientations while edge dislocations are only weakly bent. Density accumulates in the ϕ = 3π/2 direction, which is the orientation of edge dislocations required to accomodate the bending strain gradient. At the same time, the density of dislocations of other orientations decreases. The decrease is most pronounced for the ϕ = π/2 orientation, i.e. for edge dislocations of the opposite sign.

As the strain increases, we observe that the total dislocation density in the plastic regions increases strongly, as shown in Fig. 12. This increase goes along with an increase of the α 11 component of the Kröner-Nye tensor, i.e. accumulation of 'geometrically necessary' dislocations required to accomodate the increasing plastic strain gradients around the narrowing central elastic core. This is illustrated in Fig. 13.

The increase of total dislocation density that comes with the accumulation of geometrically necessary dislocations leads, according to the Taylor relation, to an increase in flow stress. As strain gradients are bigger in smaller specimens and therefore the accumulation of geometrically necessary dislocations is more pronounced, this leads to a size dependent hardening. This is illustrated in Fig. 14 which shows the scaled bending moment M/h 2 as a function of the average plastic strain γ = (1/h) +h/2 -h/2 |γ(x)| dx. While the initial flow stress is size independent, the hardening rate and the flow stress at finite strain increase with decreasing film thickness. An often-repeated argument for analyzing this type of size-dependent behavior (see e.g. [9]) runs as follows: The dislocation density in a material is the sum of a geometrically necessary density ρ G (up to a geometrical factor η equal to the strain gradient divided by the Burgers vector modulus) and a so-called 'statistically stored' density ρ S . (The latter is in fact nothing but the total dislocation density minus the geometrically necessary one, but the argument treats it like an independent physical variable.) From Taylor's relation for the flow stress one

Sfrag replacements

M = 1 • 10 -15 GPa m 2 M = 2 • 10 -15 GPa m 2 M = 5 • 10 -15 GPa m 2 -l 4 -l 2 0 + l 4 + l 2 x h = 2 µm h = 3 µm h = 5 µm h = 10 µm h → ∞ γ [-] M/h 2 [GPa]
Figure 14. Bending curves for various film thicknesses; the curve for infinite thickness has been calculated by assuming ideal plasticity and using the standard assumption of straight specimen cross-sections.

obtains

τ y = aGb √ ρ = aGb √ ρ S + ρ G = aGb ρ S + η b |∇γ| . (71)
One then makes the crucial assumption that the strain dependence of the 'statistically stored' density ρ S (γ) does not depend on the presence or absence of strain gradients. Estimating the characteristic strain gradient as γ /L where γ is a characteristic strain and L a characteristic dimension of the specimen then leads to a scaling relation for the characteristic flow stress τ :

[τ γ,L ] 2 = [τ γ,∞ ] 2 + a 2 G 2 ηb γ L , (72)
where τ γ ,∞ = aGb ρ S (γ). To test this scaling relation we have to translate the notations to our problem. The characteristic dimension L corresponds to the film thickness h, and the macroscopic flow stress corresponds to the scaled bending moment M/h 2 . As strain measure we use the average plastic strain γ as defined above. The infinite-system limit is obtained by solving the bending problem for an ideally plastic material with flow stress corresponding to the initial dislocation density. The corresponding deformation curve is denoted as [M γ /h 2 ] ∞ := lim h→∞ [M γ ,h /h 2 ]. The total dislocation density ρ in the specimen is obtained by averaging ρ(x, ϕ) over the specimen cross section and integrating over all orientations ϕ, and the 'geometrically necessary' density ρ G is obtained as the minimum density of dislocations on each of the two slip system required to accomodate the plastic strain gradient: We evaluate the component α 11 of the dislocation density tensor, average this over the specimen cross section, divide by b and multiply with the geometrical factor 1/(2 sin θ). If the above argument would apply to our simulations, a plot of [M γ ,h /h 2 ] 2 -[M γ /h 2 ] 2 ∞ versus the scaled plastic strain γ b/h should yield a straight line passing through the origin. The same should be true for a plot of the total dislocation density ρ (reduced by the infinite-system value) versus the scaled plastic strain. Our simulations, however, indicate a different behavior. Even though the 'geometrically necessary' dislocation density increases approximately linearly with scaled plastic strain, the same is not true for the total dislocation density: the increase in total density initially lags behind the increase in the 'geometrically necessary' density and only catches up at large strains (Fig. 15). As a consequence, the scaling plot of the reduced bending moment vs. the scaled plastic strain (Fig. 16) fails to produce a straight line.

What are the reasons for the discrepancy between the simulation results and the naive scaling argument? Our investigation of dislocation density evolution indicates that the linear increase of the 'geometrically necessary' density with reduced strain is not matched by a concomitant increase of the total density, and hence of the flow stress. The reason is that we allow dislocations to change their orientation, besides freely entering or leaving through the specimen surface. As a consequence, part of the increase of the 'geometrically necessary' density is not produced by additional dislocations entering through the surfaces, but rather by changes in orientation of existing dislocations and by the loss of dislocations of the 'wrong' sign. As a consequence, the total dislocation density increases, at least for small strains / large system sizes, much less than expected according to the simple scaling argument. This indicates that the assumption of a separate 'statistically stored' density that evolves in a size-independent manner is not warranted.

→ ∞ γ [-] [M γ ,h /h 2 ] 2 -[M γ /h 2 ] 2 ∞ GPa 2 < γ > b h [-]

Conclusion and Outlook

We demonstrated that ECT as a generalisation of the classical Kröner theory can capture several fundamental phenomena of the evolution of systems of curved dislocations which cannot be recovered from other dislocation based continuum theories.

In particular, the simple benchmark problem of an ensemble of expanding loops showed that the averaging process underlying ECT conserves important kinematic properties of dislocation loops. Studying the expansion of a single quasi-discrete loop served to elucidate certain numerical problems arising from the discretization of orientation space and to find appropriate remedies. We then demonstrated that the evolution of dislocation densities in orientation space may be a key factor in the description of inhomogeneous deformation processes, as changes in dislocation orientation may be an important mechanism for providing the 'geometrically necessary' dislocations needed to accomodate strain gradients. Accordingly, simulations of microbending using a simple Taylor-based flow rule yielded results that differ substantially from the predictions of gradient plasticity theories that assume strain gradients to be accomodated exclusively by the introduction of additional excess dislocations.

The present theory offers a versatile tool for treating the evolution of dislocation systems without having to resolve the discrete dynamics of all dislocation lines, or having to make arbitrary assumptions regarding the 'composition' of the dislocation arrangement. Rather, the distribution of dislocations over the various orientations follows naturally from the theory once the basic constitutive relations have been specified which relate the dislocation velocity to the dislocation pattern. In our simulations of microbending we used a most simple constitutive law, assuming local validity of the Taylor relationship but neglecting non-local contributions to the flow stress that may be associated with line curvature and/or with the piling up of dislocations of the same sign [26]. Such contributions can be easily included

Figure 1 .

 1 Figure 1. Continuous lifted curve in the configuration space. The lower, closed loop is the spatial loop, the upper line is the lifted loop with additional orientation information. The arrows attached to the lower loop indicate the spatial velocity, the arrows attached to the upper curver indicate the generalised velocity along the line. The rotational velocity ϑ is the vertical component of the generalised velocity.

2 bFigure 2 .

 22 Figure 2. Sign convention for the line orientation ϕ, tangent l and normal ν on the parameterized unit circle

2 .

 2 Discretization of the problem

( d )Figure 3 .

 d3 Figure 3. Evolution of the (spatial) total dislocation density for an expanding loop as obtained from solving Equation (50) and (51): (a) initial configuration, (b) end configuration. The loop was discretized with 60 nodes in ϕ-direction. While the total amount of density integrated over the plane is -up to a numerical error -predicted correctly (c), the loop heavily fragments. This also displays in the evolution of the numerical divergence depicted in (d).

Figure 4 .

 4 Figure 4. Discrete density distribution which corresponds to a lifted line for a bad resolution in ϕ-direction (12 nodes only) chosen deliberately for demonstration purposes. Each of these blobs is moving as a whole along the direction indicated by the arrows.

Figure 5 .

 5 Figure 5. Expansion of a continuous loop and the effect due to numerical discretisation: A lifted loop with radius r i expands by ∆r affecting the lines' inclination (the curvature) and its length. In the case of the continuouos loop expansion the line segment gets rotated (by reduction of curvature) and stretched. The numerical scheme takes care of the rotation by an evolving curvature (not shown) but otherwise can only translate the segment without stretching.

Figure 6 .

 6 Figure 6. Diagrams illustrating schematically the density evolution of a bundle of expanding dislocation loops. Graph (a) represents loops with continuous segment orientation. The other graphs refer to loops with discretized orientation: Graph (b) shows the evolution before correction, (c) after restoring line continuity by density re-distribution, and (d) with diffusive density re-distribution.

( a )Figure 7 .

 a7 Figure 7. Expansion of a circular loop with 60 nodes in ϕ-direction: projection of the density ρ (r,ϕ) into the spatial plane. The initial configuration was the same as in Fig. 3(a).

Figure 8 .

 8 Figure 8. Temporal evolution of the sum of the absolute values of div α II for different diffusive correction terms. Ideally, the sum would be zero for all steps. Due to the numerical approximation, however, it is non-zero right from the beginning and always increases. A larger amount of diffusion decreases the sum of div α II by comparison. The optimum amount of diffusion obtained for µopt is a trade-off between minimizing div α II and maintaining the correct peak density value. Relaxation of ρ yields the best result but at a very high computational cost.

5. 1 .Figure 9 .

 19 Figure 9. Microbending simulations: Investigated slip geometry. τ 0 is the bending stress (resolved shear stress in the slip systems) in the absence of plastic deformation.

Figure 10 .

 10 Figure 10. Coordinate system used in our simulations for one of the slip systems: the system is assumed to have infinite extension in y direction.

25 ρFigure 11 .

 2511 Figure 11. Dislocation density and curvature patterns in (x, ϕ) configuration space for a film thickness of h = 3µm and bending moment (per unit length in y) M = 2 • 10 -15 GPa m 2 .

2 M = 1 • 2 xM = 1 •M = 1 • 2 - 2 Figure 12 .M = 1 • 2 -M = 1 • 2 -Figure 13 .

 212112212121213 Figure 13. Accumulated plastic strain and (11) component of the Kröner-Nye tensor α for h = 3µm and different bending moments

M = 1 • 10 - 15 GPa m 2 M = 2 • 10 - 15 GPa m 2 M = 5 • 10 - 15 GPa m 2 -Figure 15 .M = 1 • 10 - 15 GPa m 2 M = 2 • 10 - 15 GPa m 2 M = 5 • 10 - 15 GPa m 2 -

 11015221015251015215110152210152510152 Figure 15. Scaling plot of the total and 'geometrically necessary' dislocation density vs reduced strain γ b/h. The upper line group is the plot for ρ , the lower line group is the plot for ρ G .

Figure 16 .

 16 Figure 16. Scaling plot of the bending curves:[M γ ,h /h 2 ] 2 -[M γ /h 2 ] 2∞ vs reduced strain γ b/h.

  (b). The growth of

	August 2, 2009 11:38 Philosophical Magazine pm-paper-revised-v1 2.0e6 2.0e6
	0.0006 0.00055	2.5e6 3.0e6 3.5e6	Philosophical Magazine 2.5e6 3.5e6 3.0e6 0.00055	15
		4.0e6		4.0e6
		4.5e6		4.5e6
	0.0008	5.0e6 6.0e6		5.0e6 6.0e6
		7.0e6		7.0e6
	0.0010		0		0
	0.0012	1000		1000
	5.0e5	2000		5.0e5 2000 1.0e6
	1.5e6	3000		1.5e6 3000 2.0e6
	2.5e6				2.5e6
	3.5e6 4.5e6			(a) Initial configuration (R = 30 nodes) F o r	3.0e6 3.5e6 4.0e6 4.5e6 5.0e6	(b) Expanded loop after 3500 steps with ∆t = 0.1 and v = 1.0 (R = 3800 nodes)
		line length	6.0e6 e P e r (c) Temporal evolution of the systems' total 7.0e6 line length e e i v number of steps numerical, ∆t = 0.1 numerical, ∆t = 0.01 analytical 1.0e6 2.0e6 3.0e6 4.0e6 5.0e6 6.0e6 number of steps 0.0004 0.0006 0.0008 0.0010 0.0012 7.0e6 mean divergence per node 0.0002 R 0 0 1000 2000 3000 0 1000 2000	3000
					w
					O n l
					y

Taylor & Francis and I.T. Consultant

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

In our case, we pragmatically consider a numerical scheme to be well-behaving within the 'working regime' if a positive half-wave can be propagated for more than 10000 steps with a velocity v = 1.0, the mesh width h = 1.0 and a time step ∆t = 0.1, such that the peak of the wave can clearly be identified and no oscillations occur.

Acknowledgement

Financial support of the European Commission under NEST-2005-PATH-COM-043386 is gratefully acknowledged.


S Civiš
email: civis@jh-inst.cas.cz

I Matulková

J Cihelka

P Kubelík

K Kawaguchi

V E Chernov

Time-resolved FTIR emission spectroscopy of Cu in the 1800-3800 cm -1 region: Transitions involving f and g states and oscillator strengths

Keywords: numbers: 52.38.Mf, 07.57.Ty, 32.30.Bv, 31.15.B- Time-resolved Fourier transform infrared spectroscopy, time decay of spectral lines, Atomic ablation, Cu, Copper, Rydberg states, oscillator strengths, line strengths, dipole transition matrix elements

    Time-resolved FTIR emission spectroscopy of Cu in the 18003800 cm1 region: transitions involving f and g states and oscillator strengths S Civiš, I Matulková, J Cihelka, P Kubelk, K Kawaguchi,

Introduction

The great advantages of Fourier transform infrared spectroscopy (FTIR), such as its constant high resolution and energy throughput, have made the IR spectral region more accessible for laboratory spectral measurements [START_REF] Nilsson | Atomic Spectroscopy and Oscillator Strengths for Astrophysical and Laboratory Plasma[END_REF]. The infrared range is becoming more and more important in astronomical research, for instance in studies of dust-obscured objects and interstellar clouds, cool objects such as discs, planets and the extended atmospheres of evolved stars, including objects at cosmological distances from the Earth (Kerber et al., 2009). Nevertheless, the powerful capacities of IR astronomy (e. g. ESOs CRyogenic Infra-Red Echelle Spectrograph, CRIRES) cannot be fully utilized without detailed spectroscopical information, first of all, on atomic line wavelengths and oscillator strengths in the IR region [START_REF] Biémont | Symposium of the International-Astronomical-Union -1st International Meeting devoted to Infrared Physics[END_REF][START_REF] Grevesse | Conference of the European-Group-for-Atomic-Spectroscopy (EGAS)[END_REF][START_REF] Pickering | International Colloquium on Atomic Spect and Oscillator Strengths (ASOS 6)[END_REF][START_REF] Jorissen | [END_REF][START_REF] Johansson | High Resolution Infrared Spectroscopy In Astronomy[END_REF]. A general drawback of the IR spectral region is the much lower number of atomic and ionic lines available (as compared to the visible and ultraviolet range) [START_REF] Ryde | [END_REF], so reporting new IR atomic lines is important.

Cu is the element (together with Zn) immediately following the iron peak, and its abundance is important in studies of scenarios of the chemical evolution of the Galaxy implied by different nucleosynthesis models (Sneden et al., 1991;Matteucci et al., 1993;[START_REF] Prochaska | [END_REF]Cunha et al., 2002;Mishenina et al., 2002;Simmerer et al., 2003;[START_REF] Bihain | [END_REF]Ecuvillon et al., 2004;[START_REF] Sobeck | First Stars III[END_REF]. Copper abundance can be analyzed from the IR spectra of cool stars [START_REF] Ryde | Atomic Spectroscopy and Oscillator Strengths for Astrophysical and Laboratory Plasma[END_REF][START_REF] Wahlgren | Cool Stars, Stellar Systems and the Sun[END_REF], metal-poor stars (Cowan et al., 2002;Honda et al., 2007) and substellar objects [START_REF] Jones | Ultralow-Mass Star Formation and Evolution Workshop[END_REF].

It should be noted that the studies of copper abundances in stars and the interstellar medium were performed mostly with optical 510.554 nm, 521.820 nm and 578.212 nm (Sneden et al., 1991;Mishenina et al., 2002;Cunha et al., 2002;Simmerer et al., 2003;Ecuvillon et al., 2004;[START_REF] Prochaska | [END_REF] or ultraviolet 327.395 nm [START_REF] Bihain | [END_REF][START_REF] Sobeck | First Stars III[END_REF] and 324.754 nm [START_REF] Sobeck | First Stars III[END_REF] Cu I lines. To our knowledge, no infrared lines of Cu I were involved in the abundance analysis except the 809.263 nm line used in a study of stellar abundances in the thick disc of the Galaxy [START_REF] Prochaska | [END_REF]. Nor have we encountered literature data about Cu I lines with wavelengths above 2 microns, and the present work is intended to supply information about Cu I infrared lines and transition moments in the 2.7-5.5 micron range. This work continues the previous studies on Au I (Civiš et al., 2010) and Ag I (Civiš et al., 2010).

Atomic copper has been attracting the great interest of spectroscopists for many decades, since its spectrum gives one of the best examples of the majority of peculiarities in atomic spectra due to a number of series perturbations observed below the first ionization threshold. Analysis of the main features of the Cu I spectrum began in the 1920s, and the first critical and exhaustive study of Cu I lines and levels was reported by Shenstone (1948). He published the spectrum from 150.4 nm to 1822.9 nm (5484-66485 cm -1 ) including both his own measurements and the previous results of other authors.

The simplest Cu I term structure originates from the closed 3d 10 core yielding 3d 10 nl j states with a total angular momentum J = j = l ± 1 2 . Other terms arise due to excitation of a 3d electron into a n ′ l ′ state. The core-excited 3d 9 nln ′ l ′ Cu I states were studied both by numerical ab initio calculations (Martin and Sugar, 1969;Carlsson, 1988) and by VUV photoabsorption [START_REF] Tondello | [END_REF]Longmire et al., 1980;Baig et al., 1992); for a more detailed review see (Baig et al., 1997). The most comprehensive reports of the closed-core 3d 10 nl j states of Cu were made by Longmire et al. (1980) The aim of the present work is the study of the Cu I spectrum in the infrared domain using FTIR spectroscopy of laser-ablated Cu plasma. The lines in this domain are due to transitions between the levels with n = 4-7; the energies of some levels (for f and g states) have not been previously reported. To identify the transitions involving these and other low-Rydberg Cu levels, we use the information on the relative intensity of the transitions obtained here by the calculation of dipole transition matrix elements (or transition dipole moments) using the Fues model potential (FMP) approach (Civiš et al., 2010). The comparison tables presented in Section 3.1 show a reasonable agreement between our calculations and the experimental and theoretical results for oscillator strengths. The dipole transition matrix elements calculated in the Section 3.1 are used in Sec. 3.2 for classification of the observed lines resulting in revised values for some terms of Cu I.

Method

The experimental setup has already been described in detail in our previous papers (Civiš et al., 2010;Kawaguchi et al., 2008). The time resolution FTIR spectra were measured using the modified Bruker IFS 120 HR spectrometer in J. Heyrovský Institute in Praha. The modification of the apparatus for the time-resolution scan of emission data was previously developed in Okayama University. The Bruker high resolution interferometer is calibrated against the internal stabilised HeNe laser with a precision around 0.001 cm -1 .

The Bruker system was equipped with an analog-digital converter (ADC 4322: Analogic, USA), which was connected to a PC containing a programmable control processor of Field Programmable Gate Array -FPGA, (ACEX 1K: Altera, USA) set up at a frequency of 33 MHz, and digital input board PCI (2172C: Interface, Japan). The data collection process and synchronization with the laser were controlled by the FPGA processor programmed by QUARTUS II 7.1, Altera. The software for data acquisition and fast FT transformation and displaying of the data were written in C ++ language.

Time-resolved FTIR spectroscopy was applied for observations of the emission arising after the irradiation of metals with a pulsed nanosecond ArF (λ = 193 nm) laser.

A high repetition rate ArF laser ExciStar S-Industrial V2.0 1000 (193 nm, laser pulse width 12 ns, frequency 1 kHz) with 15 mJ output power was focused on a rotating and linearly traversing copper target inside a vacuum chamber (average pressure 10 -2 Torr). The infrared emission (axial distance from the target 10 mm) was focused into the spectrometer using a CaF 2 (100 mm) lens. The emission was observed in the 1800-3800 cm -1 spectral region with a time profile showing maximum emission intensity at 18-20 µs after the laser shot.

For data sampling we used the so-called 1/3 sampling, where the scanner rate was set to produce a 3 kHz HeNe laser interference signal, the ArF laser oscillation was triggered, and 60 sets of time-resolved data were recorded with a preset time interval of 1 µs. Three scans were needed for a complete interferogram, and only 5 scans were coadded to improve the signal-to-noise ratio (SNR). The spectral resolution in such a procedure was about 0.1 cm -1 . We also performed measurements with only one scan but with higher resolution of about 0.017 cm -1 but worse SNR. The acquired spectra were post-zerofilled (zero filling 2, trapezoid apodization function) using the Bruker OPUS sowtware and subsequently corrected by subtracting the blackbody background spectrum. The wavenumbers, line widths and their intensities were then obtained using the OPUS internal peak picking procedure.

The Fues model potential method for the calculation of dipole matrix elements was outlined in the previous paper (Civiš et al., 2010). For Cu I the radial quantum number n r (which is equal to the number of nodes of the radial wavefunction) was assumed as:

n r = { n -4 for s, p, d-states (n ≥ 4); n -l -1 for l ≥ 3 states (n ≥ l + 1).
(1)

Results and discussion

Calculation of the dipole transition matrix elements

It is generally considered that the choice of a Coulomb approximation monoconfigurational approach (such as the FMP approach used here) is justified for transitions involving high-energy unperturbed Rydberg states, these transitions appearing mostly in IR region considered here. This approach however is reasonable for calculation of transitions involving also low-excited and even ground states. The comparison of FMP calculations with experimental and theoretical results (including ab initio calculations) for transitions in Ag was presented in the previous paper (Civiš et al., 2010). It was shown that the agreement of FMP calculations with the experiment and with the results of other calculations is satisfactory for the majority of transitions. For few transitions there were some discrepancies. However, the results of other theoretical calculations for these transitions differ between each other by some orders of magnitude and some are close to our results. From the other hand, the FMP calculations themselves are not our main aim in the present work; we use them in analyzing the relative intensities of the observed IR transitions only. To show that FMP calculation of dipole transition matrix elements are adequate for such a purpose we compare some FMP-calculated Cu oscillator strengths with the experimental and theoretical data available in the literature.

The results of such a comparison with the experiment and with other calculations are presented in Table 1. The overall agreement of our calculations with other results can be considered to be reasonable. There are some large discrepancies between our data and the values from the half-century old sources given by Allen and Asaad (1957) for 4p-ns transitions with n = 7, 8, 9. In such cases, Table 1 contains the corresponding oscillator strength values calculated in Coulomb-like approximation (Lindgård et al., 1980). We did not present here a full comparison with Lindgård et al. (1980) because of the large data reported in it. However, since the Coulomb-like approximation used by Lindgård et al. (1980) is in essence similar to FMP, the majority of the results of Lindgård et al. (1980) coincide with the FMP values with a discrepancy of not more than 20%. We demonstrate this in the Table 2 for the transitions involving the Cu I states observed in our experiment. All the uncertainties in the tables below are given in round brackets after the corresponding values and should be treated as their rightmost significant digits, e. g. 123.4(56) means 123.4 ± 5.6.

The further analysis of the observed Cu lines given in the next Section is based on FMP calculation of the line strengths S for transitions between the 3d 10 nl j states which can be easily obtained from oscillator strengths f or transition probabilities A presented in Table 2. When available, the corresponding results of (Lindgård et al., 1980) were presented for comparison; the other values have not been reported previously. The Table 2 can serve as a supplement to the tables of Lindgård et al. (1980) and to the review by Fu et al. (1995).

Lines observed

Figure 1 shows parts of the observed IR emission spectra of Cu I at 20 µs after the laser shot, when the emission intensity is maximal for almost all of the observed lines. The list of the IR lines observed for Cu I is presented in Table 3. Their full widths at half-maxima (FWHM) are calculated from fitting to the Lorentzian shape (Civiš et al., 2010;Civiš et al., 2010).

As in the previous papers (Civiš et al., 2010;Civiš et al., 2010) we have measured the emission spectrum at a different delay time, from 0 to 60 µs after the laser shot. This allows us to record the the time profiles of the observed Cu lines, i. e. their emission intensities I(t) as functions of the delay time, t. Such information can be helpful for diagnostics of the electronic state populations of neutral atoms in plumes formed by pulsed laser ablation Furusawa et al. (2004); however there are few reports of investigations on such population dynamics Rossa et al. (2009). The temporal dynamics of several lines is shown in Figure 2. While the temporal decay of some lines can be fitted, at least roughly, by an exponential function [START_REF] Hannaford | [END_REF] c Numerical Coulomb-like approximation (Lindgård et al., 1980) d Relativistic Hartree-Fock calculation with account for core polarisation effects (Migda lek and Baylis, 1978) e Relativistic Hartree-Fock calculation with model potential accounting for exchange and core polarisation (Migda lek, 1978) f Relativistic Hartree-Fock calculation with model potential and core polarisation (Migda lek and Baylis, 1979) g Calculation using the level-crossing measurement of lifetimes and configuration coupling coefficient deduced by fitting other experimental measurements of lifetimes and relative oscillator strengths (Siefart et al., 1974) h Atomic absorption measurements on flames (Lvov, 1970) i Atomic-beam absorption (Bell and Tubbs, 1970) j Critical survey of experimentally-determined oscillator strengths (Corliss, 1970) k Curves of growth in absorption measurement (Moise, 1966) l Rozhdestvenskii's hook method (Slavenas, 1966) m Atomic-beam absorption (Bell et al., 1958) n Arc emission measurements (Allen and Asaad, 1957) Table 2. FMP-calculated transition dipole moments (oscillator strengths f ik , transition probabilities A ki ) between the Cu I states observed in the present work.

I(t) = I background + I 0 exp ( - t -t 0 T ) , ( 2
)
The energies of all levels are taken from the present measurement except those for 4f and 7s levels taken from Shenstone (1948). The air wavelengths λ are calculated using these energy values.

Transition

Lower level (cm -1 )

Upper level (cm -1 )

ν (cm -1 ) λ (nm) f ik A ki (s -1 )
5f several lines display essentially non-exponential behavior including some "plateaux" or even secondary maxima at 35-50 µs after the laser shot. Their decay time, T , values are therefore estimated in Table 3 in a rough approximation; it is seen from this table that for essentially non-exponential decays the uncertainity ∆T is of the same order of magnitude as T itself. Note that this temporal decay is due to a complex combination of the collisional cascade repopulation of the emitting levels (Civiš et al., 2010) and the transfer processes in ablation products (Kawaguchi et al., 2008). The decay times T given in Table 3 is not related to the radiative lifetimes of Cu 3d 10 nl levels which are at least two orders shorter (typical orders are ∼ 10 1 -10 3 ns (Fu et al., 1995)).

After the assignment we refined the energy values for some levels involved with the classified transitions; the revised values of these energies are presented in Table 4. This procedure was similar to that used in the previous paper (Civiš et al., 2010), but unfortunately there is no high precision data on these levels, so the overall uncertainty of the refined level energy values is roughly that achieved in the old measurements (0.01-0.03 cm -1 by Shenstone (1948) and 0.07 cm -1 by Longmire et al. (1980)). Nevertheless, the difference between the energy value obtained in this work for 5f 7 2 level and that reported long ago (Shenstone, 1948) is about 2.5 cm -1 and lies beyond the uncertainty limits. Due to the better resolution of our measurements, the values listed in the Table 4 can be considered to be more accurate.

It is interesting to note that the fine-structure 5p doublet (fine-structure splitting is about 0.3 cm -1 ) is well resolved in our experiment, unlike the previous measurements (Shenstone, 1948;Longmire et al., 1980) where only a single blended line was observed. The ratio of the 5p3

Conclusion

This work continues the series of studies of emission IR spectra of metal vapors formed in ablation by a pulsed laser radiation (Civiš et al., 2010;Civiš et al., 2010) where TR FTIR spectroscopy is applied to observations of the emission arising after the irradiation of a copper target with a pulsed laser. The knowledge of Cu spectra is extended by reporting 18 lines not previously observed. All are classified as due to transitions between lowexcited states: 3d 10 ns with n = 6, 7, 3d 10 np (n = 5, 6), 3d 10 nd (n = 5, 6), 3d 10 nf (n = 4, 5, 7) and 3d 10 ng (n = 5, 6, 7).

The line classification is performed using relative line strengths expressed in terms of transition dipole matrix elements calculated with the help of the Fues model potential; these calculations show agreement with the experimental and calculated data available in the literature. In addition to these data we calculate the transition probabilities and oscillator strengths for transitions between the reported 3d 10 nl j states of Cu I.

This study reports revised values for energies of eleven known and of seven previously not reported 3d 10 nl levels of neutral Cu with n = 4-7 and l = 5-7. The newly-found 5g, 6g and 7g 2 G terms are inverted. We also record time profiles of the observed lines as functions of the delay time (0-60 µs) after the laser shot with maxima of emission intensity at 18-20 µs after the shot. Some lines display singleexponential temporal decay, while the decay of other lines is essentially non-exponential and demonstrates some "plateaux" (or even secondary maxima) at 40-50 µs after the laser shot. The approximate decay time for different lines varies in the 5-15 µs range. Longmire et al. (1980)

  (for ns states up to n = 41 and np states up to n = 57) and by MacAdam et al. (2009) (for s, p, d, f and g states for n = 23-28).

Figure 1 .Figure 2 .

 12 Figure 1. Observed emission IR spectrum of Cu I with a detailed structure of the most prominent lines around 2480-2530 cm -1 . Both the inset and the main spectrum correspond to the resolution of 0.017 cm -1 .

  Figure 3. The 5p1 2 , 3 2 -6s 1 2 doublet lines of Cu I in two different resolutions. The 5p 1 2 -6s 1 2 component (with center of gravity at 3465.481 cm -1 ) is split due to hyperfine structure.

Table 1 .

 1 Comparison of FMP-calculated oscillator strengths for transitions between some 3d 10 nl j states of Cu I with the values from other sources

	Transition FMP	Other values
		(this work)
	4s 1 2 -4p	0.223	0.22175 a ; 0.215(6) b ; 0.2631 c ; 0.214 d ; 0.224 e ; 0.215 f ; 0.220(15) g ; 0.22 h ; 0.2230 j ; 0.323(4) l ; 0.153 k ; 0.16(4) m ; 0.0158 n
	4s 1 2 -4p	0.449	0.44756 a ; 0.439(12) b ; 0.5296 c ; 0.432 d ; 0.453 e ; 0.434 f ; 0.432(28) g ; 0.43(2) i ; 0.4460 j ; 0.66 l ; 0.322 k ; 0.31(3) m ; 0.0316 n
	4p 1 2 -5s 4p 3 2 -5s 4p 1 2 -6s 4p 3 2 -6s 4p 3 2 -6d 5 2 4p 1 2 -7s 4p 1 2 -8s 4p 1 2 -9s 4p 3 2 -9s 4p 1 2 -4d 3 2 4p 3 2 -4d 5 2 4p 3 2 -4d 3 2 4p 1 2 -5d 3 2 4p 3 2 -5d 5 2 4p 3 2 -5d 3 2	0.171 0.172 0.0134 0.0131 0.0465 0.00419 0.00191 0.00105 0.00102 0.592 0.538 0.0598 0.131 0.118 0.0131	0.151 e ; 0.1690 j 0.153 e ; 0.1730 j 0.0136 e ; 0.00909 j ; 0.0148 n 0.0136 e ; 0.0131 j ; 0.0157 n 0.04441 c ; 0.00378 n 0.004064 c ; 0.00723 n 0.001899 c ; 0.00245 n 0.001053 c ; 0.000251 n 0.001027 c ; 0.00056 n 0.51106 a ; 0.551 d ; 0.546 e ; 0.4880 j 0.46727 a ; 0.503 d ; 0.499 e ; 0.4650 j 0.05201 a ; 0.056 d ; 0.0555 e ; 0.0614 j 0.1238 c ; 0.117 e ; 0.0932 j 0.1114 c ; 0.107 e ; 0.0791 j 0.01236 c ; 0.0119 e

a Solution of one-electron Schrödinger equation with Hartree-Slater potential corrected by core polarisation and spin-orbital interaction

(Curtis and Theodosiou, 1989)

b Laser-induced fluorescence from sputtered metal vapors

Table 3 .

 3 Experimental Cu I lines and their identification. The decay time, T , was calculated by exponential fitting of the measured time profiles of the corresponding lines. Time profile demonstrates significant deviation from the exponential decay b The decay curve has essentially non-exponential form with a plateau or secondary maxima; τ value is absent or roughly approximate

	Wavenumber	Intensity	SNR FWHM	Decay	Identification
	(cm -1 )	(arb.		(cm -1 )	time
		units)			(µs)
	1887.307(6) 1935.313(2) 2163.890(16) 2171.118(18) 2176.426(16) 2179.011(3) 2494.8098(3) 2497.7750(3) 2513.814(4) 2517.4511(3) 2521.0550(3) 2865.233(2) 3110.955(4) 3465.481(4) 3465.8044(7) 3837.402(12) 3840.376(15)	1.30×10 4 6.58×10 4 3.42×10 2 1.64×10 2 2.16×10 2 7.68×10 4 2.98×10 5 3.77×10 5 4.65×10 3 1.03×10 5 1.58×10 5 2.32×10 4 2.06×10 4 8.46×10 4 1.70×10 5 7.49×10 3 9.64×10 3	17 63 6.2 4.6 5.6 58 63 121 8.4 78 129 37 27 10 23 9.5 11	0.083(16) 0.106(5) 0.104(46) 0.074(55) 0.077(43) 0.065(10) 0.038(1) 0.041(1) 0.048(13) 0.051(1) 0.049(1) 0.062(5) 0.087(11) 0.089(13) 0.048(2) 0.193(36) 0.235(47)	12.1(44) b 5.88(221) b 16.2(41) a 9.50(266) a 12.8(34) b 9.06(479) a 11.7(42) b 11.3(36) b 11.5(41) b 15.0(57) b 14.9(57) b 7.93(167) b 8.53(214) b 5.06(87) b 5.25(90) b 8.51(209) a 14.7(190) a	6p 3 2 -7s 6s 1 2 -6p 5f 7 2 -7g 5f 5 2 -7g 6d 5 2 -7f 6s 1 2 -6p 4f 5 2 -5g 4f 7 2 -5g 5d 5 2 -5f 5d 3 2 -5f 5d 5 2 -5f 6p 1 2 -6d 6p 3 2 -6d 5p 1 2 -6s 5p 3 2 -6s 4f 5 2 -6g 4f 7 2 -6g

a

Table 4 .

 4 Revised values of some levels of Cu I

	Term	Energy (cm -1 ) Other sources
	3d 10 5p 1 2 3d 10 5p 3 2 3d 10 5d 3 2 3d 10 5d 5 2 3d 10 5f 5 2 3d 10 5f 7 2 3d 10 5g 7 2 3d 10 5g 9 2 3d 10 6s 1 2 3d 10 6p 1 2 3d 10 6p 3 2 3d 10 6d 3 2 3d 10 6d 5 2 3d 10 6g 7 2 3d 10 6g 9 2 3d 10 7f 7 2 3d 10 7g 7 2 3d 10 7g 9 2	49383.263(23) 49383.26 a 49382.949(14) 49382.95 a 55387.621(11) 55387.668 a 55390.569(9) 55391.292 a 57905.041(14) 57905.2 a , 57905.23 b 57911.090(12) 57908.7 a 57924.610(30) This work 57924.075(30) This work 52848.752(9) 52848.749 a 55027.763(26) 55027.74 a , 55027.713 b 54784.081(21) 54784.06 a , 54784.073 b 57893.028(24) 57893.05 a 57895.084(24) 57895.1 a 59267.202(33) This work 59266.676(34) This work 60071.510(30) This work 60076.159(23) This work 60074.980(20) This work

a

Shenstone (1948) b

Acknowledgments

This work was financially supported by the Grant Agency of the Academy of Sciences of the Czech Republic (Grants No. IAA400400705 and KAN 100500652).