Particle-based approach to the Eulerian distortion field and its dynamics

The Eulerian distortion field is an essential ingredient for the continuum modeling of finite elastic and inelastic deformations of materials; however, its relation to finer levels of description has not yet been established. This paper provides a definition of the Eulerian distortion field in terms of the arrangement of the constituent microscopic particles, which is beneficial for fundamental studies as well as for the analysis of computer simulations, e.g., molecular dynamics simulations. Using coarse graining and nonequilibrium thermodynamics, the dynamics of the Eulerian distortion field is examined in detail and related to the underlying dynamics of the particles. First, the usual kinematics of the distortion and the known expression for the Cauchy stress tensor are recovered. And second, it is found that the Mandel stress and the plastic deformation-rate tensor in the natural configuration constitute the relevant force–flux pair for the relaxation of the distortion. Finally, the procedure is illustrated on two examples, namely on an amorphous solid and on a crystalline solid with one slip system.

Inelastic behavior is intrinsic to practically all materials. While the reasons for such behavior differ significantly between materials, there are also common aspects, particularly when it comes to the general structure of the modeling approaches. In principle, the current stress state depends on the entire history of imposed deformation (e.g., by way of a memory function/kernel), which is at the basis of Boltzmann's principle of superposition [3,4] used to describe linear viscoelasticity. However, if nonlinearity comes into play, a different approach is often taken. A prominent alternative is to decompose the total deformation into its elastic and inelastic contributions. In the field of solids, this is realized by the multiplicative decomposition of the total deformation gradient F, according to Kröner [5] and Lee and Liu [6,7], which is a cornerstone for modeling finite deformation plasticity (see, e.g., [8][9][10][11][12][13]): It uses F = F e · F p , where F e and F p represent the elastic and plastic parts of the total deformation gradient, respectively. The quantity F p maps line elements from the reference configuration to an intermediate stress-free configuration, from which F e maps these line elements to the current configuration. Conceptually, F e can be determined by the instantaneous (to avoid relaxation effects) unloading of a deformed sample. The main goal of such approaches is the formulation of an evolution equation for F e that contains both an affine contribution (representative of the kinematics of F) and a contribution that describes the rate of plastic deformation (representative of the inelastic dynamics of F p ); notably, the plastic deformation F p itself is not involved, as long as, for example, strain/work hardening is not considered. A discussion of some critical aspects of the multiplicative decomposition, particularly in relation to the intermediate configuration, can be found in [14]. A concept related to that of the intermediate configuration is the so-called natural reference state [15]; for a discussion on similarities and differences between the two, the reader is referred to [15]. In either case, the elastic deformation F e is the fundamental dynamic variable for modeling the material behavior under finite deformation.
Instead of the elastic deformation F e , one could equivalently use the so-called distortion A, which is the inverse of F e , A = F −1 e . This is the approach taken in [16][17][18]. Although the distortion and deformation gradient are equivalent in the purely elastic case, where A = F −1 , distortion becomes advantageous in the case of dissipation, where the distortion only reflects the elastic part of the deformation gradient, A = F −1 e . Dissipation can actually break the continuum compatibility conditions: the configuration from which the elastic part of the deformation gradient is mapping, which is called the natural configuration in the distortion-related literature, ceases to exist as a global configuration [19][20][21][22]. This happens, for instance, in the case of a screw dislocation [23], where the natural configuration has a discontinuity, and thus, we cannot take derivatives with respect to it, and F e does not exist globally. In the case of dislocations, the multiplicative splitting can be substantiated mathematically by projection of the deformation tensors to a regular part and a part that contains the Dirac distribution [14]. In contrast to the natural configuration, the current configuration always exists, so we can take derivatives with respect to it, and A exists. If the distortion is curl-free, its contour integrals give the natural configuration, that represents the equilibrium state upon instantaneous unloading [19,20]. If, however, the distortion is not curl-free (for instance in case of dislocations, where the curl of the distortion equals the dislocation density tensor [21,[24][25][26][27]), then the result of the integration depends on the particular path chosen, which means that the natural configuration is defined only locally.
Corresponding to the dynamics of the elastic deformation (e.g., see [10,11]), the dynamics of the distortion can be written in the general form [17] D t A = −A · (∂v/∂r) where D t denotes the material derivative, v is the velocity field, and p L and pĽ stand for the plastic deformationrate tensors in the natural and current configurations, respectively. If one disregards inelastic deformation, it is obvious that by virtue of its definition, the deformation gradient is in principle a Lagrangian field, with the reference position as dummy field variable [28,29]. In contrast, the distortion is quite naturally an Eulerian field, with the current position as dummy field variable. This being said, it is of course possible to, once defined, re-interpret the deformation gradient as well as its elastic and plastic parts as Eulerian fields (e.g., see [30,31]). Nevertheless, when aiming at an Eulerian formulation for modeling finite deformation inelastic behavior, the Eulerian distortion simplifies the description because it involves only one mapping (from the current to the reference configuration); in contrast, for instance, the Eulerian deformation gradient is the composition of two mappings (mapping from the reference to the current configuration, but position evaluation in the current configuration).
In this paper, the goal is to express the continuum Eulerian field of distortion in terms of the arrangement of the constituent microscopic particles of the material and to derive the distortion dynamics Eq. (1) in relation to the underlying microscopic-particle dynamics. While such an approach has already been taken earlier in relation to the deformation gradient [31], it is still lacking for the distortion field, to the best of our knowledge. Linking a particle-based description to continuum quantities is a prototypical example of coarse graining and multiscale modeling [32,33]. Doing so is not only beneficial for a fundamental understanding of the continuum model, it also helps to devise efficient hierarchical modeling strategies.
The developments in this paper concern the deformation in D-dimensional Euclidean space, e.g., D = 3 for bulk materials or D = 2 for plane, flat domains. It is pointed out for completeness that one is not dealing here with curved spaces, e.g., a curved surface such as the surface of a soap bubble, because there the dynamics involves contributions that are perpendicular to the tangent spaces to the surface [34].
Before we proceed, a word on notation: Discrete particles are enumerated by Greek indices (α, β,...), while coordinates are labeled with Latin indices: lowercase for the current configuration (i, j,...) and uppercase for the natural configuration (I , J ,...). For example, coordinates of particle α in the current and natural configurations are given by r i α and R I α , respectively, while A I j denotes a component of the distortion tensor. In all that follows, the notation of tensor calculus according to Ricci-Curbastro and Levi-Civita [35,36] is used. This implies that upper and lower indices are distinguished corresponding to contravariant and covariant components of tensors. Throughout this paper, the upper/lower positioning of indices is used, and it is always understood-without explicit mention-that one can transition between the two representations by multiplication with the appropriate metric. Summations over particle indices are always spelled out, while for summation over coordinates the Einstein summation convention is used, and summations are only permitted over index pairs of which one index is upper and the other is lower. For practical convenience, for the developments and derivations in this paper, Cartesian coordinates with Euclidean metric-given by the Kronecker delta (where the distinction between upper and lower indices is irrelevant)-are used, without loss of generality. However, specifically in Sect. 6 the main results of the paper are reformulated for generalized (curvilinear) coordinates. Finally, for functions f of macroscopic position, the short-hand notation f ≡ f (r) and f ≡ f (r ) is often used for better readability.
The paper is organized as follows: Sect. 2 starts with deriving a best-fit expression for the distortion field in terms of discrete particle positions on a finer level, which is then employed to examine the reversible (affine) and irreversible (plastic) dynamics of the distortion by direct calculation. In Sect. 3, the General Equation for the Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) framework is introduced, with special emphasis on its systematic procedure for statistical-mechanics-based coarse graining. This procedure is then employed to examine in more detail the reversible dynamics (Sect. 4) and the irreversible dynamics (Sect. 5) of the distortion, respectively. In Sect. 6, the main results of this paper are presented in covariant form, for generalized coordinates. Conclusions are drawn in Sect. 7. The main achievements are summarized as "Result" in the respective sections.
2 Particle-based approach to distortion 2.1 Particle-based definition of the distortion An essential ingredient in relating particle-based approaches, e.g., molecular dynamics simulations, to continuum mechanics is to provide a particle-based definition of the macroscopic measure of deformation. For an introduction and discussion of best-fit (and similar) approaches for obtaining microscopic definitions of strain tensors, see [37]. A frequently used procedure is that developed by Falk and Langer [38], and very similar ones have been used by other researchers (e.g., see [39][40][41]), where the configurations of particles in two different states are related by a best-fit deformation gradient. Since only particles are considered which are in the neighborhood of a certain particle α, this best-fit deformation gradient is in that sense local (from a macroscopic perspective) and in particular it is a best-fit deformation gradient for that specific particle α. In order to use this for the definition of a macroscopic deformation gradient field, one can add up all these contributions, where each contribution gets a weight δ(r α − r), in complete analogy to how one would define the mass density field for given particle masses m α . This is the approach taken earlier [31]; however, it has the unwanted side effect that the thereby defined macroscopic deformation gradient has a density-like character, and an additional step for proper normalization is needed (see [31] for details).
Here, we take another approach for finding a particle-based expression for the distortion. Rather than defining the distortion for the neighborhood for a certain particle and then processing this further into a field quantity, we right-away select only those particles that are in the vicinity of a certain macroscopic position r. Specifically, we propose that the distortion at macroscopic position r is given by the minimization of the error, in adaptation of the procedure of Falk and Langer [38], where w α = w(r α − r) with a localization function w and both particle index summations (α and β) run over all particles. The positions of the particles in the current configuration are denoted by r α , whereas their positions obtained after instantaneous unloading, i.e., in the natural configuration, are given by R α ; their connector vectors are denoted by r αβ = r β − r α and R αβ = R β − R α , respectively. The introduction of the weights, w α w β , ensures that only particle pairs are considered for which both of the particles are in the current configuration close to the macroscopic position r. By V I 2 , we mean the square of the norm of the vector V, V I 2 = V I V I , where V I and V I denote the covariant and contravariant components of the vector, respectively. It is noted that, for the definition of the distortion, it is not required that the localization function w be normalized; it is only later-in relation to the definitions of the densities of mass and momentum (Sect. 4)-that normalization to unity is needed. Throughout this paper, it is assumed that the normalization function satisfies w(−x) = w(x).
In terms of mechanics, the natural configuration is stress-free and obtained by unloading the current configuration infinitely rapidly (to avoid any relaxation effects, which would be nonelastic) to a state at ambient isotropic pressure. In terms of molecular simulations, unloading infinitely rapidly implies that this process should not be performed according to physical dynamics according to, e.g., a molecular dynamics (MD) simulation. Rather, one needs to consider the landscape of the proper energy function, in this case the Gibbs free energy with pressure as a variable, and use steepest descent in phase space to arrive at a local minimum in this landscape. The so-obtained configuration is what we understand as the natural configuration on the level of the discrete constituent particles. It is closely related to what is called an "inherent structure" of the landscape of the potential energy [42,43]. A specific procedure for the examination of local minima in the landscape of the Gibbs free energy has been discussed in detail in [44,45].
Minimization of the error Eq. (2) with respect to (w.r.t.) the components of the distortion leads to the best-fit distortion, given by with Not only based on its definition, but also from the expressions in Eq. (3), it is evident that the right side of the distortion refers to the current configuration, whereas the left side refers to the natural configuration.
The following aspects of this expression for the distortion are noteworthy: (i) This expression for the distortion does not depend on the normalization of the localization function w, as anticipated. (ii) More importantly, this expression for the distortion is automatically normalized properly: It does not suffer from having density-like character as in [31], while nevertheless being local, i.e., it takes into account only those particles that are in the vicinity of the macroscopic position r. (iii) This expression for the distortion becomes ill-defined only if there are no (or rather: not sufficiently many) particles in the vicinity of r, which is, however, reasonable and in agreement with the macroscopic scenario.

Result 1 Equation
(3) is a particle-based definition of the distortion in an Eulerian setting (i) with control over the size of the zone-of-influence (by way of the localization function w), (ii) that is normalized properly, and (iii) that is ill-defined only in case of a lack of sufficient particles in the zone-of-influence.

Reversible dynamics of the distortion: direct calculation
The reversible dynamics of the distortion describes the kinematics in the current configuration, while the natural configuration (i.e., the unloaded state) is not affected. Therefore, the reversible dynamics of the distortion Eq. (3) has its origin exclusively in the evolution of the particle positions in the current configuration, r α . To proceed, the particle dynamics is assumed to be well approximated bẏ with v being the velocity field. This approximation should be precise enough due to the presence of the kernel w(r α − r) in Q and P, the product of which is the distortion A. The reversible evolution of the distortion is then given byȦ I k =Q I j P jk + Q I jṖ jk , where the dot refers to the rate-of-change due to the particle dynamics described by Eq. (4). Using the expressions forQ I j andṖ jk derived in "Appendix A," the reversible evolution of the distortion thus becomeṡ The third contribution on the right-hand side (r.h.s.) is negligible compared to the other ones, because both particles α and β must be close to each other (due to the kernels w α and w β ) and because, in this region of interest, the distortion has been obtained by a best-fit procedure, which implies Therefore, the reversible evolution of the distortion simplifies tȯ which agrees with [22,46] and is fully compatible with the kinematics of the total deformation gradient (e.g., see [29]). In other words, the particle-based definition of distortion Eq. (3) leads to the correct reversible evolution of the distortion. This result is summarized in the following:

Result 2
The particle-based definition of the distortion, Eq. (3), the approximate particle evolution, Eq. (4), and the local equilibrium between the particle positions in the current and natural configurations, Eq. (6), lead to the standard reversible contribution to the evolution of the distortion, Eq. (7).

Irreversible dynamics of the distortion: preliminary comments
As can be seen above, the reversible dynamics affects distortion on its right side, i.e., on the side that relates to the current configuration. In contrast, the irreversible dynamics is expected to affect the distortion on its left side, i.e., on the side that relates to the natural configuration. This expectation is in line with Eq. (3), since the left side of the distortion is related to the particle positions in the natural configuration, {R α }, via Eq. (3b). The notion that the irreversible (i.e., inelastic) dynamics is expressed most clearly in the natural configuration is not new, but has been put forward explicitly in [11,15]. To make this more concrete, one can proceed as follows. Based on Eq. (4) for the reversible particle dynamics in the current configuration, one findsṙ i αβ = (∂v i /∂r j )r j αβ . By analogy, one may thus approximate the purely irreversible dynamics of the particle connector vectors in the natural configuration bẏ with p L being a local plastic deformation-rate tensor in the natural configuration. If one assumes, for illustration purposes only, that p L is the same for all particle pairs (α, β), it can be shown by straightforward calculation that the irreversible contribution to the dynamics of the distortion Eq. (3) is equal to the second term on the r.h.s. of Eq. (1a). While this merely serves as a first motivation, the effect of the particle dynamics in the natural configuration on the irreversible evolution of the distortion will be studied more thoroughly in Sect. 5.

Method: nonequilibrium thermodynamics and statistical mechanics
To approach the dynamics of the distortion from a wider perspective, we make use of nonequilibrium thermodynamics. While a multitude of modeling frameworks have been developed in this field over the years, we specifically choose the General Equation for the Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) framework [47][48][49]; the main reason for this choice is that this framework also comes with a procedure for coarse graining, i.e., for establishing links between models on different levels of description, which is the main focus of this paper in relation to the distortion. In the following, only a very concise summary of the framework is outlined as needed for this paper; for more details, the reader is referred to the original papers on GENERIC [47,48] and its coarse graining procedure [50], the comprehensive books [49,51], or applications of this framework to the case of fluid and solid mechanics, e.g., [30,31,52]. Let X denote the set of fundamental variables, for which a closed set of evolution equations is sought; with reference to the specific case of interest in this paper, X could contain the momentum density, the mass density, the temperature (or density of internal energy or entropy), and the distortion. According to GENERIC, the evolution equations are given by [47][48][49] with total energy E, entropy S, Poisson operator L, and friction operator M. The two contributions on the r.h.s. of Eq. (9) are reversible (e.g., affine deformation, elasticity) and irreversible (e.g., plasticity) in nature, respectively. The dot (·) does not only represent summation over discrete indices in the sense of different elements of X , but may also stand for integration over continuous labels in the case of field theories, e.g., as what will be discussed in the remainder of this paper. The two operators must satisfy various properties: On the one hand, (L1) L must be anti-symmetric, (L2) L must satisfy the Jacobi identity (to ensure time-structure invariance), and (L3) the derivative of entropy, δS/δX , must be in the null-space of L, i.e., L · (δS/δX ) = 0; on the other hand, (M1) M must be Onsager-Casimir symmetric, (M2) M must be positive semi-definite, and (M3) the derivative of energy, δ E/δX , must be in the null-space of M, i.e., M · (δ E/δX ) = 0. It can be shown that these six conditions imply the conservation of energy E and that the entropy S is a non-decreasing function of time.
The fluctuation-dissipation theorem [53][54][55] states that, whenever there is dissipation, there are in principle also fluctuations present; whether the fluctuations are practically relevant is case specific. In GENERIC, the dissipation originates from the presence of the friction operator M. Incorporation of (Wiener process based) fluctuations in a thermodynamically consistent manner extends the GENERIC from Eq. (9) to the following stochastic differential equation (SDE) [47,49,50], using the Itô interpretation of stochastic calculus [56,57], with Boltzmann constant k B and where dW stands for multicomponent white noise, more precisely, for the increment of a multicomponent Wiener process [56,57]. The "strength" B of the fluctuating contribution in Eq. (10) is related to the friction operator by way of which ensures that the fluctuation-dissipation theorem is respected. The deterministic case Eq. (9) can be recovered from the fluctuating one, Eq. (10), in the limit k B → 0, while leaving the building blocks E, S, L, and M unchanged. As already mentioned, the GENERIC framework can be used also for establishing a relation between different levels of description. Imagine to coarse grain from a fine level with variables X 1 to a coarser level with variables X 2 . The instantaneous value of the level 2 quantities, 2 , can be expressed in terms of level 1 quantities, 2 = 2 (X 1 ), and the coarse-grained variables are eventually given by X 2 = 2 , where . . . denotes the average with an appropriate distribution function [49,50]. In analogy to conventional statistical mechanics, one can derive expressions for the energy and entropy, E 2 and S 2 , on level 2 based on level 1 information; since this is not pursued in this paper, it is not presented in this concise summary, and we refer to [49,50] for details on this specific aspect.
The reversible and irreversible dynamics on level 1 is also represented on level 2. Concretely, this is reflected by the relations [31,49,50] where the contractions (·) work between the level 1 operator and X 1 . Furthermore, there might be an additional dissipative contribution on level 2, since some features of the dynamics that were explicitly resolved on level 1 may be seen on the coarser (and slower) level 2 as rapid, uncontrollable noise. In this case, there is a second contribution to the friction operator on level 2 [31,49,50], which involves the correlation of the fluctuating contribution to the dynamics of 2 ,˙ f 2 , at different instances in time. The quantity τ denotes the separating time-scale between levels 1 and 2. If the correlation function decays sufficiently fast, Eq. (14a) can be written in the form [49] where τ f 2 stands for the fluctuation-related increment of 2 in the time interval τ . The prefactor on the left-hand side (l.h.s.) of Eq. (14b) is a short-hand notation for the following: Since X 2 in general contains several variables, M 2 has a matrix structure. The prefactor on the l.h.s. in Eq. (14b) means that each element of M 2 is to be multiplied by a factor that depends on the parity under time reversal (ε = ±1) of the two respective variables concerned, which is directly related to the Onsager-Casimir reciprocal relations [58][59][60]. In particular, if the two variables concerned have the same parity, the prefactor equals 2, which is the case for what will be discussed in this paper.
A slightly different formulation of GENERIC, where the friction operator is replaced by a dissipation potential [47], is advocated in [61]. Although the formulation with a friction operator is more general than the formulation with dissipation potential [62] (for instance, because the operator may not always be symmetric), the dissipation potential is appealing because of an intimate connection with the principle of large deviations [63,64]. It remains unclear, however, which of the formulations should be preferred in general. Here, we choose the formulation with the friction operator because it contains a standardized procedure of statisticalmechanics-type coarse graining [49,50].
The following two sections focus on the derivation of both the reversible and irreversible parts of the dynamics of the distortion from the level of discrete microscopic particles.

Full set of variables and generating functionals
Before proceeding with the derivation of the distortion dynamics based on coarse graining (specifically deriving the operators L and M), the general form of the energy E and entropy S is briefly discussed. The set of variables is chosen to consist of the momentum density m, the mass density ρ, the entropy density s, and the distortion A, i.e., the classical hydrodynamic fields, augmented by the distortion. All densities are of Eulerian type, i.e., w.r.t. the current configuration. Instead of the entropy density, another thermal variable could have been chosen, e.g., the density of internal energy e, or the absolute temperature T ; however, the entropy density s turns out to be particularly useful, specifically when formulating the reversible dynamics.
The corresponding expressions for E and S are then written as where dV r denotes the volume element in the current configuration. The concrete expression for the density of internal energy e is not specified here, but it could be derived by systematic coarse graining if so desired. There are two assumptions implicit in Eq. (16): (i) Next to the explicit kinetic-energy contribution, the remaining (internal) energy does not depend on the momentum density field; (ii) the density of internal energy depends on the fundamental variables only through their local values, i.e., non-local effects are neglected. This choice is made for convenience only, rather than for necessity. Based on Eq. (16), the functional derivatives of the energy and entropy are found to be which generate the dynamics, according to Eqs. (9) and (10). In passing, it is mentioned that the velocity (being the time derivative of position) has an upper index, and conversely, its conjugate variable, the momentum (density), has a lower index.

Derivation of Poisson operator by coarse graining
In order to derive the Poisson operator L 2 by way of coarse graining, Eq. (12), microscopic expressions are needed for the instantaneous values of the level 2 quantities, i.e., 2 , in terms of level 1 quantities, X 1 ; for the latter, we choose the positions and momenta of all particles. While the expression for the distortion has already been provided above, Eq. (3), the expressions for the mass density and momentum density are given by with mass m α and momentum p α of particle α, w α = w(r α − r), and where the localization function w must be normalized to unity, wdV r = 1. A microscopic expression for the entropy density is not specified at this point; the entropy density will be discussed further below. The Poisson bracket on level 1, which is related directly to the Poisson operator L 1 , is given by the well-known expression for classical mechanics [65], Correspondingly, the Poisson operator on level 2, Eq. (12), which will be denoted by L instead of L 2 from now on for simplicity, can be written as [31,49,50] where and ♦ are two variables on level 2, i.e., from Eq. (15).
It is obvious from Eq. (20) that the corresponding element of the Poisson operator L is nonzero only if at least one of the two quantities and ♦ depends on the particle momenta. In view of the microscopic expressions for the momentum density m(r) (Eq. (18b)), the mass density ρ(r) (Eq. (18a)), and the distortion A(r) (Eq. (3)), it follows that only the m-related column and row of L have nonzero entries. For these entries, explicit calculations based on Eq. (20) lead to the following expressions: where the approximations in Eqs. (21a) and (21b) consist in using that the localization function w can be approximated by a Dirac delta function on level 2, i.e., on the macroscopic scale. The element which couples distortion with momentum, Eq. (21c), is the most difficult one to calculate, and therefore, only the result is shown here, while the derivation is found in "Appendix B." It is noted that the expressions for L m i (r),m j (r ) and L ρ(r),m j (r ) agree with what has been used earlier in a purely continuum formulation [48,49] and by coarse graining [31]. What remains to be discussed are the elements of the Poisson operator related to the entropy density, s(r). What is the reversible dynamics of such a field? There are two ways to find the dynamics. The first way relies on the requirement that the overall entropy cannot be changed by the reversible dynamics, which results in a degeneracy condition for the Poisson operator. The other way uses coarse graining similar to what we have done for the other fields; however, this would require a definition of the entropy density in terms of some microscopic variables. In particular, one could start with the Liouville N -particle distribution function and formulate the local entropy density in terms of the one-and two-particle distribution functions within the BBGKY hierarchy [66,67]. While this latter approach is more in the spirit of coarse graining as compared to the first one, it would make the derivation of the Poisson operator more complicated and thus distract from the main goal of this paper, which is a particle-based approach to the distortion. Therefore, we shall follow the first (simpler) procedure. The conservation of entropy can be formulated in the same way as the conservation of mass, which means that the component of the Poisson operator coupling the entropy and momentum is [49] in complete analogy to Eq. (21b) for the mass density. In summary, the Poisson operator L is completely specified by the elements in Eq. (21) and using its anti-symmetry property, whereas all other elements are zero. It can be shown that this Poisson operator agrees with the Poisson bracket used earlier for distortion dynamics [46,68]. Therefore, the derivation of the Poisson operator is summarized in the following result: with ∂ t ≡ ∂/∂t and ∂ i ≡ ∂/∂r i , and v i = δ E/δm i . The Cauchy stress tensor is given by with hydrostatic pressure The evolution of the distortion, Eq. (22d), has been written in two equivalent forms: The first one shows explicitly that the evolution depends on the curl of the distortion, while the second highlights the split into advection and deformation. Although the distortion can be thought of as the gradient of the mapping from the current configuration to the natural configuration, it does not need to be curl-free because the natural configuration does not need to exist as a globally defined manifold [16,20,22], see also Sect. 1. Moreover, if we dropped the curl terms, we would violate Galilean invariance of the evolution equation [46]. The curl terms are also important for the validity of the Jacobi identity, which would cease to be unconditionally valid without them [46]. In order to demonstrate that the stress tensor Eq. (23) is symmetric, it is used that the thermodynamic potential e must not change with rotations in the current configuration. This implies that e depends on the distortion only by way of the inverse of the (elastic) right Cauchy-Green deformation tensor (C e ), with contravariant components h I K ≡ C −1,I K e = A I j A K j . With this, the Cauchy stress tensor Eq. (23) becomes which makes the symmetry manifest also for the second contribution. The evolution of the distortion as described by Eq. (22d) agrees with earlier results in the literature [22,46,68] and is also fully compatible with the kinematics of the deformation gradient [29,52]. Furthermore, the stress tensor Eq. (23) coincides with earlier distortion-based formulations [22,46], and after rewriting the stress tensor in terms of the (elastic part of the) deformation gradient, Eq. (23) is found to be in agreement with the results in [11,[28][29][30][31]52].

Derivation of friction operator by coarse graining
If the fine level of description, i.e., that of the discrete particles, has no irreversible dynamics on its own, the friction operator M on the coarse level consists only of a contribution of the fluctuation-form Eq. (14), whereas a contribution of the form Eq. (13) is absent. As far as the fluctuations are concerned, it is noted that, due to the conservation of total energy, the fluctuations must satisfy the constraint [49,50] δ In the following, the focus is exclusively on the relaxation of the distortion as the only irreversible contribution to the dynamics, thereby neglecting the conduction of heat and other irreversible processes. It is reasonable to assume that the relaxation of the distortion does not affect the balances of mass and momentum. Therefore, Eq. (25) implies If one introduces the abbreviation T = ∂e/∂s, and if one assumes that the condition Eq. (26) is fulfilled pointwise, one can express the fluctuations of the entropy in terms of fluctuations of the distortion, Using the general fluctuation formula for the friction operator, Eq. (14a), together with relation Eq. (27) between the fluctuations, all nonzero elements of the coarse-grained friction operator can be related to its By way of its construction, it follows that the friction operator satisfying the restrictions in Eq. (28) complies with the degeneracy condition M · (δ E/δX ) = 0 (with δ E/δX evaluated at r and performing an integration over r ), i.e., that the total energy is not affected by the irreversible dynamics, as required. In order to derive by coarse graining the remaining element of the friction operator, M A I j A K l , one needs some information about the irreversible dynamics on the particle level. In the present case, this is related to the purely irreversible dynamics of {R α } in the natural configuration, in which the plastic effects are represented most naturally (see also [11,15]). The way forward depends in general on the form of this dynamics, in particular whether only small or also large fluctuations play a role. Since the natural configuration is obtained after instantaneously unloading the current configuration (under certain boundary conditions), {R α } denotes a local minimum in the energy landscape [69]. The evolution of the natural configuration can therefore be imagined as a jump process between local minima [70,71]. If these jumps are rare events, the coarse graining procedure needs to reflect that aspect properly [72,73]. In the following, it is assumed for the purpose of illustration that the dynamics of {R α } takes the form of a SDE with multiplicative noise, using the Itô interpretation of stochastic calculus [56,57], where both a α and b αβ can be functions of the particle arrangement. The symbol dW β stands for increments of Wiener processes, with the properties Equation (30b) for the second moment is the component-wise representation of dW α (t)dW β (t ) = dtdt δ(t − t )δ αβ 1.
The fluctuation contribution to the distortion dynamics can be obtained by way of τ is given by the second term on the r.h.s. of Eq. (29). With this, the coarse graining expression for the friction operator, Eq. (14b), becomes where we have made use of τ W K α τ W β,L = τ δ αβ δ K L , by virtue of Eq. (30b). Equation (31) describes how the particle dynamics in the natural configuration, Eq. (29), gives rise to dissipation on the coarse-grained level. It is by way of taking the average that Eq. (31) will turn into a closed-form expression for M A I j A K l in terms of the coarse-grained variables Eq. (15).
For the sake of turning Eq. (31) into a closed-form expression by analytical calculation, a special case is considered, namely Writing the coefficients c I J K in terms of information in the current configuration is chosen for convenience during the further steps of the analysis; intuitively, one can think of the dependence on r αβ as a dependence on R αβ by way of Eq. (6), and additionally, the dependence on r α allows for incorporation of a dependence on the distortion evaluated at the position of particle α, e.g., to encode stress-dependent kinetics. About Eq. (32), it is noted that it depends on particles α and β only, which is physically intuitive. However, when performing the double summations over particles in Eq. (31), the presence of the localization function w enables the replacement of particle positions by macroscopic positions r and r , respectively. In this way, it can be shown that the expression in the first bracket on the r.h.s. of Eq. (31) can be approximated by −c I M J (r γ − r, r) Q J m , and similarly, the second bracket in Eq. (31) can be approximated by −c K M L (r γ − r , r ) Q Ln . With this and using Eq. (3) for the distortion, Eq. (31) simplifies significantly. If the statistical average is taken w.r.t. a distribution that is microcanonical in the distortion (i.e., proportional to δ( A J j − A J j )), A J j inside of the average may be replaced by A J j and taken out of the average, which leads to with the correlation function defined by In the following, two examples for the correlation function will be examined. Example 1 For illustration purposes, the most simple ansatz is that there is a complete lack of correlation at different positions and in different directions, where the scalar function z accounts for the overall strength of the correlation. Two possible realizations of where in both cases the vector field d must satisfy Some possible choices for d are discussed in "Appendix C". Example 2 As a second example, crystal plasticity is considered; for an introduction, see, e.g., [74]. The plastic slip system is described in terms of the glide direction s and the glide-plane normal n, both having unit length. Both of these vectors are physically relevant and should therefore appear in the expression for the third-rank tensor c I J K ; a third vector is needed, which is called q for now. Since it is wanted that the vectors s and n contribute to the index structure of the correlation functions Eq. (34) in Eq. (33) and eventually to the coarse-grained dynamics, this latter vector q will have to be placed in between s and n. In other words, the indices of q are contracted in the correlation function Eq. (34). A class of propositions for the correlation function is thus given by where q is a yet unknown vector and s has been placed on the left in order to allow R-fluctuations only in the glide direction s. If one chooses, e.g., q = d defined via Eq. (37), the ansatz Eq. (38) leads to where D denotes the number of space dimensions. The consequences that these two examples for the correlation function have for the irreversible dynamics will be discussed in the following.

Irreversible contributions to the evolution equations
According to the GENERIC Eq. (9), the irreversible contribution to the evolution equations given by the entropy derivative in Eq. (17), evaluated at r , and the friction operator given by the elements in Eq. (28) with Eqs. (33) and (34), with an integration over r , are where the driving force for plastic flow and the plastic deformation-rate tensor in the natural configuration (see also Eq. (1a)) are given by It can be shown that Eq. (41) is actually equal to the Mandel stress [75], apart from an isotropic (hydrostatic) contribution; however, that difference can be traced back to our use of ρ and A as independent variables. In Eq. (40c), the entropy inequality ∂ t s ≥ 0 is obviously satisfied, due to the symmetry and positivity of the integrand by virtue of the properties of Z I J K L according to Eq. (34). The general structure of the irreversible contributions to the evolution equations, Eq. (40), is compatible with earlier work about dissipative distortion dynamics [22,46].
The significance of the Mandel stress [75] as the relevant quantity in relation to large-deformation plasticity modeling is well documented in the literature [9,[11][12][13]15,31,76]. The Mandel stress and the plastic deformation-rate tensor in the natural configuration form the relevant force-flux pair in the context of nonequilibrium thermodynamics, a fact that has been noted and exploited earlier [9,[11][12][13]15,31,76] and which is also represented by the second equation in Eq. (40c) and by Eq. (42).
The concrete form of the irreversible dynamics of the distortion depends on the example considered, as shown in the following.
where z ≥ 0 may be a function of the Mandel stress. It may be useful to rewrite the expression Eq. (43) in terms of the plastic deformation-rate tensor in the current configuration, which is often used in the modeling of amorphous solids. Based on the equivalence of the different representations of plastic deformation in Eq. (1), one finds Here, c ik = A I i A I k are the covariant components of the (elastic) Cauchy deformation tensor, the latter being also known as the Piola tensor or Finger tensor; furthermore, it can be shown that c ik = B −1 e,ik where B e is the elastic left Cauchy-Green deformation tensor. Equation (44) shows that only if the Cauchy stress tensor commutes with the Cauchy deformation tensor, pĽi j is directly proportional to the (non-hydrostatic part of the) Cauchy stress tensor in terms of its index structure. This latter case is standard when modeling the plasticity of isotropic materials, see, e.g., [22,30,77,78].
Example 2: Using the correlation function Eq. (39), the plastic deformation-rate tensor in the natural configuration becomes where the expression in parenthesis is the resolved shear stress [8,74]. This form for the plastic deformationrate tensor is in agreement with standard literature [8,10,74,79], where the kinetic coefficient z ≥ 0 may be a function of the Mandel stress.

Incorporation of incompressibility of irreversible deformation
In the above discussion about plasticity, no special attention has been paid to the evolution of the determinant of the distortion A. While the mass density ρ has been chosen as an independent variable and remains unaffected by the irreversible dynamics in Eq. (40), it may nevertheless be desirable to discuss the conditions under which also the determinant of the distortion remains constant. This aspect is considered in detail in the following. For simplicity, the reversible contribution to the dynamics is neglected in this subsection.
To start, the case is considered where the dynamics of the distortion with increment d A I j is deterministic, i.e., where fluctuations are absent. Requiring that the determinant of the distortion remains constant can be cast into the condition 0 = d (ln(det A)). With ∂ (ln(det A)) /∂ A I j = A −1, j I , this condition can be written in the form For illustration purposes, consider the irreversible dynamics of the distortion represented in Eq. (1), expressed in terms of the plastic deformation-rate tensor, either in the natural or current configuration. The condition Eq. (46a) implies that the trace of both plastic deformation-rate tensors must vanish. For the further procedure, it is useful to express the condition Eq. (46a) for the case where the dynamics is represented in the GENERIC Eq. (9); for this case, it can be shown that Eq. (46a) assumes the form It is mentioned in passing that another procedure for the discussion of incompressibility has been introduced in terms of projection operators [30,31]; however, the conditions in Eq. (46) are more useful for the discussion here.
If the dynamics of the distortion, d A I j , does contain fluctuations, the above discussion needs to be amended. In general, it can be stated that whenever there is a friction operator, there can in principle also be fluctuations, as shown in Eq. (10). In the present context, the fluctuations in the distortion dynamics are related to the fluctuations of the particle arrangement in the natural configuration. In line with the treatment in Sec. 5.1, it is thus assumed that the distortion dynamics can be described in terms of an SDE based on Wiener processes [56,57]. Starting from the deterministic dynamics of the distortion, the extension toward the SDE with fluctuations is completely analogous to the treatment in [80], and the reader is referred to that paper for the details. Once the SDE for the distortion dynamics is determined, the change in ln(det A) can be obtained as follows. Using Itô's formula of stochastic calculus [57,81], one finds where ∂(ln(det A))/∂ A I j = A −1, j I and ∂ 2 (ln(det A))/∂ A I j ∂ A K l = −A −1, j K A −1,l I have been used; in the second contribution on the r.h.s. of Eq. (47), only the fluctuations need to be included and reduced according to Itô's scheme (see Table 3.1 in [57]). The crucial observation-after a straightforward though lengthy calculation-is that Eq. (46b) is not only the necessary condition for incompressibility in the deterministic case, but also it is the necessary and sufficient condition to ensure incompressibility on average, even in the presence of fluctuations. One may thus make use of the condition Eq. (46b) as a means of verifying or enforcing the incompressibility in the dynamics of the distortion. The question is now how the incompressibility condition on the friction operator, Eq. (46b), is to be implemented on the level of the particle dynamics in the natural configuration. To that end, it is instructive to consider Eq. (33) for the distortion dynamics. Multiplication of Eq. (33) with A −1, j I A −1,l N (multiplication with A −1,l N is just for convenience and without loss of generality, since A is invertible) leads to the necessary conditions (note the implicit summations), where Eq. (49b) follows directly from Eq. (49a) and the symmetry of the correlation function, Eq. (34). The conditions Eq. (49) have direct consequences for the relation Eq. (42) between the plastic deformationrate tensor in the natural configuration and the Mandel stress: First, Eq. (49a) results in tr( p L) = 0, which is the well-known formulation to ensure the incompressibility of plastic deformation (see, e.g., [11,15]). And second, Eq. (49b) implies that p L is insensitive to isotropic contributions in the Mandel stress, which in turn means that the true driving force for incompressible plastic deformation is only the deviatoric part of the Mandel stress (see, e.g., [11,15]). This section is completed by examining the effects of incompressibility on the two examples introduced earlier.
Example 1: The correlation function specified in Eq. (35) does not satisfy the incompressibility condition Eq. (49a). In order to encode incompressibility in the correlation function, Eq. (35) must be amended as follows: as can be verified by direct calculation. What remains to be determined are concrete expressions for c I J K that satisfy Eq. (50). To that end, we make use of the vectors d introduced earlier, satisfying the condition Eq. (37). As can be shown by direct calculation, the following two expressions are solutions to Eq. (50), as a modification of Eq. (36): For both of these expressions, c I J I = 0, and therefore, the condition Eq. (49a) is fulfilled in a strong sense, i.e., even without the summation over all particles, and without averaging.
Example 2: According to Eq. (38), it follows directly since the glide direction s is always perpendicular to the glide-plane normal n. In turn, this implies that the condition Eq. (49a) is automatically fulfilled, and the corresponding friction operator element M A I j A K l thus does account properly for the incompressibility of the plastic deformation due to slip, as expected. Point 1: Wherever Kronecker deltas appear with both indices up or both indices down, respectively, they need to be replaced by the appropriate metric: δ i j → g i j , δ i j → g i j , δ I J → G I J , and δ I J → G I J , where the so-called Eulerian metric g and the so-called material metric G denote the metrics in the current and natural configuration, respectively. Furthermore, these metrics are also to be used whenever indices are raised or lowered in the respective configuration. While the presence of the Eulerian metric helps to preserve the covariance of the equations, the material metric is necessary for instance to calculate the energy. The two (Eulerian and material) metrics act on different manifolds, which means that they are different in general [28]. If the manifolds have neither curvature nor torsion, using the Levi-Civita connection [82], we can choose such coordinate systems that the metrics are constant (by Frobenius theorem [83]) and equal to the unit matrices [84]. The material metric has already been introduced in [85] for modeling of pre-stressed materials and in [86] for space-time compatibility reasons. Point 2: Partial derivatives w.r.t. a spatial coordinate, ∂ i ≡ ∂/∂r i , need to be expressed as covariant derivatives, involving the Christoffel symbols, see [35,36,87] in general and also [11] in the context of mechanics. Obviously, the covariant derivative reduces to the partial derivative in the case of Cartesian coordinates. Point 3: Even if working in generalized coordinates, it is recommended to interpret densities w.r.t. the unit volume element in the Cartesian setting. However, this implies that, when integrating the density for obtaining the corresponding extensive quantity, the (square root of the determinant of the) metric needs to be taken into account.
Point 4: It may be sufficient to keep the definition of the distortion Eq. (3) in the Cartesian formulation, for the following reason. Due to the localization function, the region of interest is small on a macroscopic scale; in this region of interest, a local Cartesian coordinate system can be introduced as a reasonable approximation.
It is noted that particularly for Point 1 and Point 2 it proves beneficial that systematic use has been made of the notation of tensor calculus according to Ricci-Curbastro and Levi-Civita [35,36] throughout this paper and that Kronecker deltas have been used with appropriate positioning of indices.

Conclusion
In this paper, we have introduced a particle-based definition of the Eulerian distortion field by using a bestfit procedure in a certain zone-of-influence, where the latter can be controlled by a localization function (see Sect. 2). Specifically, the best-fit distortion Eq. (3) is defined as the minimizer of the error Eq. (2) that measures the discrepancy between distances in the natural and current configurations, respectively.
While a simplified analysis of the dynamics of the distortion is pursued in Sect. 2 by direct calculation, there is a need for a more thorough analysis in the form of systematic coarse graining from the microscopic dynamics, in particular to shed more light on the irreversible dynamics of the distortion. In other words, the simplified point of view of Sect. 2 has been made more precise by using the GENERIC framework and its procedure for coarse graining, as summarized in Sect. 3. Within GENERIC, evolution equations are composed of a reversible (Hamiltonian) part and an irreversible (dissipative) part, where all relevant building blocks can be expressed in terms of liner-level information. In this way, the evolution equations for the distortion and the densities of mass, momentum, and entropy have been derived by coarse graining in Sect. 4 and Sect. 5. With regard to the reversible contribution to the evolution equations (Sect. 4), the underlying Poisson operator Eq. (21) has been derived from the canonical Poisson operator for classical particles.
In Sect. 5, the emergence of the irreversible dynamics of the distortion has been examined in detail. Specifically, the relaxation of the distortion-representative of inelastic/plastic deformation-has been linked to the underlying dynamics of the constituent particles, Eq. (29), which is expressed by the friction operator Eq. (28) with Eqs. (33) and (34). Irrespective of the material considered, it has been established that the Mandel stress is the relevant driving force for the relaxation of the distortion. The coarse graining procedure has been illustrated w.r.t. the irreversible dynamics with the help of two examples, an amorphous isotropic solid and a crystal with one slip system. Finally, it has been examined how the incompressibility of distortion relaxation is expressed not only on the coarse-grained continuum level, but also on the fine level of the constituent particles.
In the following, four perspectives for future work are identified and discussed briefly. Outlook 1: The results of this paper can be used not only for analyzing the results of particle-based simulations (e.g., molecular dynamics simulations of atoms or of slightly coarse-grained, so-called united, atoms) in terms of extracting the field of distortion as such, but also for studying the rapid contributions to the fine-scale dynamics for eventually obtaining the friction operator (see Sect. 5 for details). In this respect, it is crucial to examine whether the fluctuations can be represented adequately in terms of stochastic differential equations (SDE) as discussed in Sect. 5, or whether there are rather rare but large events that prevail [88], where the latter would require other means of coarse graining (e.g., see [72,73]). An equally important question is the one of time-scale separation between the different levels of description; a lack of time-scale separation (as observed, e.g., in dislocation systems [89]) complicates the coarse graining significantly.
Outlook 2: In order to account for substantial inhomogeneities of the deformation, it has been proposed to let the strain energy depend not only on a certain measure of deformation itself, but also on its gradient, see, e.g., [90]. If the measure of deformation is given by the Eulerian distortion, this implies that the strain energy depends not only on the distortion A I j , but also on its partial derivative w.r.t. the spatial coordinates, In principle, H I jk could be calculated as the partial derivative of the distortion A I j , once the latter has been determined according to Sect. 2.1. However, particularly when discrete-to-continuum relations are of interest [90], it may be preferable to obtain H I jk directly from the arrangement of the particles, in analogy to Eq. (3) for the distortion A I j . This can actually be achieved by following the procedures developed in [41] or in [91,92], as detailed in "Appendix D." While following [41] leads to an expression for H I jk that is inherently symmetric, H I k j = H I jk , the approach taken in [91,92] is more general, see "Appendix D" for details. In either case, it might be interesting to extend the results for the dynamics of A I j derived in this paper to the case where both A I j and H I jk (as derived in "Appendix D") are independent dynamic variables. Outlook 3: It is mentioned in Sect. 1 that in the presence of irreversible deformation, the distortion is in general not compatible, i.e., not curl-free. In the case of dislocations, the curl of the distortion is related to the dislocation density tensor [21,[24][25][26][27] and the Burgers tensor [16,46]. Such a relation has actually been employed successfully in [91,92] for detecting dislocations based on configurations of the discrete particles, i.e., atoms. One may wonder whether even for a wider class of systems, the curl of the distortion is a useful tool for detecting carriers of plastic deformation. With Eq. (53), the curl of the distortion (expressed in Cartesian coordinates) reads with Levi-Civita's permutation symbol ε i jk (taking values −1, 0, 1); the extension to generalized coordinates requires simply to introduce the prefactor 1/ √ g on the r.h.s. of Eq. (54) [36,93], where g denotes the determinant of the metric in the current configuration. In line with Outlook 2, there are two possibilities for the application of Eq. (54) on the level of the particles: Either the spatial derivatives of the distortion are calculated after the distortion A I j has been obtained by the best-fit procedure in Sect. 2.1, or one makes use of H I jk that itself is obtained directly by way of a best-fit procedure (as described in "Appendix D"). In the latter case, the approach following [41] results in a vanishing curl (due to H I k j = H I jk ) and is therefore insufficient for the discussion of incompatibility. Instead, the approach following [91,92] should be used. However, while [91,92] propose a two-step (i.e., sequential) procedure for the calculation of the best-fit measure of deformation and its derivative, it is worth pursuing options for a one-step (i.e., simultaneous) procedure, as mentioned briefly in the discussion section of [92]. Whatever procedure is used for the determination of the (non-vanishing) curl of the distortion, it is emphasized that no a priori assumptions about the defect structure enter that determination; rather, it processes directly the arrangements of particles in the current and natural configurations. If the so-obtained mapping between these two configurations turns out to be incompatible, one may use the curl of the distortion as a tool for detecting defects, particularly dislocations, along the lines of [91,92].
Outlook 4: As mentioned in Sect. 1, this paper deals with flat spaces. In contrast, when studying spaces with curvature, additional care needs to be taken and amendments to the procedure developed in this paper may be required. This is because one cannot find a parametrization such that the components of the metric are equal to the Kronecker delta everywhere.
where we have used ∂w α /∂r i α = −∂w α /∂r i . For the other contribution, we note thaṫ P jk (r) = −P jl P −1,lm P mk .
For the calculation of (P −1,lm ) , one can show that, using similar steps as above and again with Eq. (4), Making use of ∂ m l (r ) /∂ p γ,k = δ k l w γ and noting that A I j does not depend on the particle momenta, Eq. (20) can be used to calculate the Poisson operator element L A I j (r),m l (r ) , based on the three contributions in Eq. (B4). The first contribution to L A I j (r),m l (r ) is obtained by inserting Eq. (B4b) in Eq. (20). That intermediate result can be processed further by the following steps: (i) We approximate w γ = w(r γ − r ) by (ii) it can be shown that the contribution related to w(r − r ) vanishes for symmetry reasons (α ↔ γ ); (iii) the remaining expression (with prefactor 2) can be rewritten once with original and once with interchanged summation indices (α, γ ), the combination of which can be simplified (iv) with the aid of Eq. (3c) and then with Eq. (3a). Eventually, one obtains On a macroscopic scale, w(r − r ) can be approximated by the Dirac delta-function δ(r − r ). Using A I l = A I l , one thus finds The second contribution to L A I j (r),m l (r ) is obtained by inserting Eq. (B4c) in Eq. (20), which leads to where we have used ∂w γ /∂r l γ = −∂w γ /∂r l . It can be shown by straightforward calculation, by using the multiplicative decomposition Eq. (3a) and a rule analogous to Eq. (A2), that the expression in the angular brackets is identical to ∂ A I j /∂r l , i.e., Again replacing w(r − r ) by δ(r − r ) on the macroscopic scale and with A I j = A I j , one obtains The third contribution to L A I j (r),m l (r ) is obtained by inserting Eq. (B4d) in Eq. (20), which leads to Following the same argument as for the second contribution above, one can approximate w(r γ −r ) ≈ w(r−r ), which leads to Since the summation runs over all values of α and γ , each index pair occurs twice with identical weight, w α w γ ; however, for these two occurrences the corresponding expressions in the parenthesis have opposite sign and thus cancel each other out. Therefore, In summary, the expression for the Poisson operator element that couples the distortion and the momentum density is given by Clearly, if the inhomogeneity of the deformation is significant on small scales, Eq. (D18) is to be preferred over Eq. (2). The goal is to minimize the error Eq. (D18a) where the tensor U −1 is the inverse of the tensor U given by i.e., with U ilmn U −1 mn jk = δ i j δ l k . Finally, Eq. (D21) can be inserted in Eq. (D20) to obtain the final expression also for the distortion.

D.2 Procedure 2
In [91,92], another procedure has been proposed for the determination of the derivatives of the deformation gradient. When applying this procedure to the case of the Eulerian distortion in the following, [91,92] are followed closely; only the main steps are presented.
After having determined the microscopic expression for the distortion A I j as described in Sect. 2.1, its spatial derivative can be determined as follows: The difference in the distortion evaluated at the positions of two particles, could also be expressed in terms of H I jk r k αβ , if the two particles are sufficiently close to each other. While these two expressions for the difference in distortion are in general not identical for all particle pairs for one and the same H I jk , one can determine H I jk as the minimizer of the error and P nk given by Eq. (3c).