text:
# References

Statistical Multifragmentation in Thermodynamical Limit:
An Exact Solution for Phase Transitions

K.A. Bugaev1,2, M.I. Gorenstein1,2, I.N. Mishustin1,3,4 and

W. Greiner1

1 Institut für Theoretische Physik, Universität Frankfurt, Germany

2 Bogolyubov Institute for Theoretical Physics, Kyiv, Ukraine

3 Kurchatov Institute, Russian Research Center, Moscow, Russia

4 Niels Bohr Institute, University of Copenhagen, Denmark

## Abstract

An exact analytical solution of the statistical multifragmentation model is found in thermodynamic limit. Excluded volume effects are taken into account in the thermodynamically self-consistent way. The model exhibits a 1-st order phase transition of the liquid-gas type. An extension of the model including the Fisher’s term is also studied. The possibility of the second order phase transition at or above the critical point is discussed. The mixed phase region of the phase diagram, where the gas of nuclear fragments coexists with the infinite liquid condensate, is unambiguously identified. The peculiar thermodynamic properties of the model near the boundary between the mixed phase and the pure gaseous phase are studied. The results for the caloric curve and specific heat are presented and a physical picture of the nuclear liquid-gas phase transition is clarified.

Key words: Nuclear matter, 1-st order liquid-gas phase transition,

mixed phase thermodynamics

PACS: 21.65.+f, 24.10. Pa, 25.70. Pq

Nuclear multifragmentation is one of the most interesting and widely discussed phenomena in intermediate energy nuclear reactions. The statistical multifragmentation model (SMM) (see  and references therein) was recently applied to study the liquid-gas phase transition in nuclear matter . Numerical calculations within the canonical ensemble exhibited many intriguing peculiarities of the finite multifragment systems. However, the investigation of the system’s behavior in the thermodynamic limit was still missing. Therefore, there was no rigorous proof of the phase transition existence, and the phase diagram structure of the SMM remained unclear. Previous numerical studies for the finite nuclear systems (the canonical and microcanonical ensembles) led to unjustified (and sometimes wrong) statements concerning the nuclear liquid-gas phase transition in the thermodynamic limit. In our recent paper  an exact analytical solution of the SMM was found within the grand canonical ensemble which naturally allowed to study the thermodynamic limit. The self-consistent treatment of the excluded volume effects was an important part of this study. In this letter we investigate the peculiar thermodynamic properties near the boundary between the mixed phase and the pure gaseous phase. New results for the caloric curve and the specific heat are presented and a physical picture of the nuclear liquid-gas phase transition in SMM is clarified. This physical picture differs from the one advocated in the previous numerical studies.

In the SMM the states of the system are specified by the multiplicity sets $`\{n_k\}`$ ($`n_k=0,1,2,\mathrm{}`$) of $`k`$-nucleon fragments. The partition function of a single fragment with $`k`$ nucleons is : $`\omega _k=V\left(mTk/2\pi \right)^{3/2}z_k`$, where $`k=1,2,\mathrm{},A`$ ($`A`$ is the total number of nucleons in the system). $`V`$ and $`T`$ are, respectively, the volume and the temperature of the system, $`m`$ is the nucleon mass. The first two factors in $`\omega _k`$ originate from the non-relativistic thermal motion and the last factor, $`z_k`$, represents the intrinsic partition function of the $`k`$-fragment. For $`k=1`$ (nucleon) we take $`z_1=4`$ (4 internal spin-isospin states) and for fragments with $`k>1`$ we use the expression motivated by the liquid drop model (see details in Ref. ): $`z_k=\mathrm{exp}(f_k/T),`$ with the fragment free energy

$$f_k=[W_\mathrm{o}+T^2/ϵ_\mathrm{o}]k+\sigma (T)k^{2/3}+\tau T\mathrm{ln}k.$$
(1)
Here $`W_\mathrm{o}=16`$ MeV is the bulk binding energy per nucleon, $`T^2/ϵ_\mathrm{o}`$ is the contribution of the excited states taken in the Fermi-gas approximation ($`ϵ_\mathrm{o}=16`$ MeV) and $`\sigma (T)`$ is the temperature dependent surface tension which is parameterized in the following form:

$$\sigma (T)=\sigma _\mathrm{o}[(T_c^2T^2)/(T_c^2+T^2)]^{5/4},$$
(2)
with $`\sigma _\mathrm{o}=18`$ MeV and $`T_c=18`$ MeV ($`\sigma =0`$ at $`TT_c`$). The last Fisher’s term in Eq. (1) with dimensionless parameter $`\tau `$ is introduced for generality. The canonical partition function (CPF) of the ensemble of nuclear fragments has the following form:

$$Z_A^{id}(V,T)=\underset{\{n_k\}}{}\underset{k=1}{\overset{A}{}}\frac{\omega _k^{n_k}}{n_k!}\delta (A\underset{k}{}kn_k).$$
(3)
The model defined by Eqs.(1,3) with $`\tau =0`$ was studied numerically in Refs. . This is a simplified version of the SMM since the symmetry-energy and Coulomb contributions are neglected. However, its investigation appears to be very important for understanding the physics of multifragmentation.

In Eq. (3) the nuclear fragments are treated as point-like objects. However, these fragments have non-zero proper volumes and they should not overlap in the coordinate space. In the Van der Waals excluded volume approximation this is achieved by replacing the total volume $`V`$ in Eq. (3) by the free (available) volume $`V_fVb_kkn_k`$, where $`b=1/\rho _\mathrm{o}`$ ($`\rho _\mathrm{o}=0.16`$ fm-3 is the normal nuclear density). Therefore, the corrected CPF becomes: $`Z_A(V,T)=Z_A^{id}(VbA,T)`$.

The calculation of $`Z_A(V,T)`$ is difficult because of the constraint $`_kkn_k=A`$. This difficulty can be partly avoided by calculating the grand canonical partition function:

$$𝒵(V,T,\mu )\underset{A=0}{\overset{\mathrm{}}{}}\mathrm{exp}\left(\mu A/T\right)Z_A(V,T)\mathrm{\Theta }(VbA),$$
(4)
where the chemical potential $`\mu `$ is introduced. The calculation of $`𝒵`$ is still rather difficult. The summation over the sets $`\{n_k\}`$ in $`Z_A`$ cannot be performed analytically because of the additional $`A`$-dependence in the free volume $`V_f`$ and the restriction $`V_f>0`$. The problem can be solved by introducing the so-called isobaric partition function (IPF) which is calculated in a straightforward way (see details in Refs. ):

$$\widehat{𝒵}(s,T,\mu )_0^{\mathrm{}}𝑑V\mathrm{exp}(sV)𝒵(V,T,\mu )=\frac{1}{s(s,T,\mu )},$$
(5)
where

$`(s,T,\mu )`$ $`=`$ $`\left({\displaystyle \frac{mT}{2\pi }}\right)^{3/2}[z_1\mathrm{exp}\left({\displaystyle \frac{\mu sbT}{T}}\right)`$ (6)
$`+`$ $`{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}k^{3/2\tau }\mathrm{exp}\left({\displaystyle \frac{(\nu sbT)k\sigma k^{2/3}}{T}}\right)],`$
with $`\nu \mu +W_\mathrm{o}+T^2/ϵ_\mathrm{o}`$. In the thermodynamic limit $`V\mathrm{}`$ the pressure of the system is defined by the farthest-right singularity, $`s^{}(T,\mu )`$, of the IPF $`\widehat{𝒵}(s,T,\mu )`$

$$p(T,\mu )T\underset{V\mathrm{}}{lim}\frac{\mathrm{ln}𝒵(V,T,\mu )}{V}=Ts^{}(T,\mu ).$$
(7)
The study of the system’s behavior in the thermodynamic limit is therefore reduced to the investigation of the singularities of $`\widehat{𝒵}`$.

The IPF (5) has two types of singularities: 1) the simple pole singularity defined by the following equation $`s_g(T,\mu )=(s_g,T,\mu )`$; 2) the singularity of the function $``$ itself at the point $`s_l(T,\mu )=\nu /Tb`$ where the coefficient in linear over $`k`$ terms of the exponent in Eq. (6) is equal to zero.

The simple pole singularity corresponds to the gaseous phase where pressure $`p_gTs_g`$ is determined by the following transcendental equation: $`p_g(T,\mu )=T(p_g/T,T,\mu )`$. The singularity $`s_l(T,\mu )`$ of the function $``$ defines the liquid pressure: $`p_l(T,\mu )Ts_l(T,\mu )=\nu /b.`$ Here the liquid is represented by an infinite fragment (condensate) with $`k=\mathrm{}`$.

In the region of the $`(T,\mu )`$-plane where $`\nu <bp_g(T,\mu )`$ the gaseous phase is realized ($`p_g>p_l`$), while the liquid phase dominates at $`\nu >bp_g(T,\mu )`$. The liquid-gas phase transition occurs when the two singularities coincide, i.e. $`s_g(T,\mu )=s_l(T,\mu )`$. As $``$ in Eq. (6) is a monotonously decreasing function of $`s`$ the necessary condition for the phase transition is that the function $``$ is finite in its singular point $`s_l`$. At $`\tau =0`$ this condition requires $`\sigma (T)>0`$ and, therefore, $`T<T_c`$. Otherwise, $`(s_l,T,\mu )=\mathrm{}`$ and the system is always in the gaseous phase as $`s_g>s_l`$. As one can see from Eq.(6) the convergence properties of $`(s,T,\mu )`$ depend significantly on the Fisher’s exponent $`\tau `$ in the vicinity of the critical point where the surface term vanishes. In what follows we mainly concentrate on the case $`\tau =0`$. Other possibilities which appear at $`\tau >0`$ are discussed shortly. Their detail study can be found in Ref. . Here we only note that Eqs. (5, 6) represent an exact solution of the Fisher’s droplet model  where additionally the effects of excluded volume are incorporated.

The baryonic density $`\rho `$ in the liquid and gaseous phases is given by the following formulae, respectively:

$$\rho _l\left(p_l/\mu \right)_T=1/b,\rho _g\left(p_g/\mu \right)_T=\rho _{id}/(1+b\rho _{id}),$$
(8)
where the function $`\rho _{id}`$ is the density of point-like nuclear fragments with shifted, $`\mu \mu bp_g`$, chemical potential:

$`\rho _{id}(T,\mu )`$ $`=`$ $`\left({\displaystyle \frac{mT}{2\pi }}\right)^{3/2}[z_1\mathrm{exp}\left({\displaystyle \frac{\mu bp_g}{T}}\right)`$ (9)
$`+`$ $`{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}k^{5/2}\mathrm{exp}\left({\displaystyle \frac{(\nu bp_g)k\sigma k^{2/3}}{T}}\right)].`$
At $`T<T_c`$ the system undergoes a 1-st order phase transition across the line $`\mu ^{}=\mu ^{}(T)`$ defined by the condition of coinciding singularities: $`s_l=s_g,`$ i.e., $`p_l=p_g`$. The phase transition line $`\mu ^{}(T)`$ in the $`(T,\mu )`$-plane corresponds to the mixed liquid and gas states. This line is transformed into the finite mixed-phase region in the $`(T,\rho )`$-plane shown in Fig. 1. The baryonic density in the mixed phase is a superposition of the liquid and gas baryonic densities: $`\rho =\lambda \rho _l+(1\lambda )\rho _g,`$ where $`\lambda `$ ($`0<\lambda <1`$) is the fraction of the system’s volume occupied by the liquid inside the mixed phase. Similar linear combinations are also valid for the entropy density $`s`$ and the energy density $`\epsilon `$ with $`(i=l,g)`$ $`s_i=\left(p_i/T\right)_\mu `$ , $`\epsilon _i=T\left(p_i/T\right)_\mu +\mu \left(p_i/\mu \right)_Tp_i`$.

Inside the mixed phase at constant density $`\rho `$ the parameter $`\lambda `$ has a specific temperature dependence shown in Fig. 2: from an approximately constant value $`\rho /\rho _\mathrm{o}`$ at small $`T`$ the function $`\lambda (T)`$ drops to zero in a narrow vicinity of the boundary separating the mixed phase and the pure gaseous phase. This corresponds to a fast change of the configurations from the state which is dominated by one infinite liquid fragment to the gaseous multifragment configurations. It happens inside the mixed phase without discontinuities in the thermodynamical functions.

An abrupt decrease of $`\lambda (T)`$ near this boundary causes a strong increase of the energy density as a function of temperature. This is evident from Fig. 3 which shows the caloric curves at different baryonic densities. One can clearly see a leveling of temperature at energies per nucleon between 10 and 20 MeV. As a consequence this leads to a sharp peak in the specific heat per nucleon at constant density, $`c_\rho (T)(\epsilon /T)_\rho /\rho `$, presented in Fig. 4. A finite discontinuity of $`c_\rho (T)`$ arises at the boundary between the mixed phase and the gaseous phase. This finite discontinuity is caused by the fact that $`\lambda (T)=0`$, but $`(\lambda /T)_\rho 0`$ at this boundary (see Fig. 2).

The negative values of the specific heat shown in Fig. 4 appear due to the parameterization of the surface tension (2): its second derivative with respect to temperature generates a negative contribution in the vicinity of $`T_c`$ for all densities. This feature of the model is unphysical and it has to be modified.

It should be emphasized that the energy density is continuous at the boundary of the mixed phase and the gaseous phase, hence the sharpness of the peak in $`c_\rho `$ is entirely due to the strong temperature dependence of $`\lambda (T)`$ near this boundary. Moreover, at any $`\rho <\rho _\mathrm{o}`$ the maximum value of $`c_\rho `$ remains finite and the peak width in $`c_\rho (T)`$ is nonzero in the thermodynamic limit considered in our study. This is in contradiction with the expectation of Refs.  that an infinite peak of zero width will appear in $`c_\rho (T)`$ in this limit. Also a comment about the so-called “boiling point” is appropriate here. This is a discontinuity in the energy density $`\epsilon (T)`$ or, equivalently, a plateau in the temperature as a function of the excitation energy. Our analysis shows that this type of behavior indeed happens at constant pressure, but not at constant density! This is similar to the usual picture of a liquid-gas phase transition. In Refs.  a rapid increase of the energy density as a function of temperature at fixed $`\rho `$ near the boundary of the mixed and gaseous phases (see Fig. 3) was misinterpreted as a manifestation of the “boiling point”.

The results presented in Figs. 1-4 are obtained for $`\tau =0`$. New possibilities appear at non-zero values of the parameter $`\tau `$. At $`0<\tau 5/2`$ the qualitative picture remains the same as discussed above, although there are some quantitative changes. For $`\tau >5/2`$ the condition $`(s,T,\mu )<\mathrm{}`$ is also satisfied in the singularity point $`s_l(T,\mu )`$ for all $`T>T_c`$ where $`\sigma (T)=0`$. Therefore, the liquid-gas phase transition extends now to all temperatures. Its properties are, however, different for $`\tau >7/2`$ and for $`\tau <7/2`$ (see Fig. 5). If $`\tau >7/2`$ the gas density is always lower than $`1/b`$ as $`\rho _{id}`$ is finite. Therefore, the liquid-gas transition at $`T>T_c`$ remains the 1-st order phase transition with discontinuities of baryonic density, entropy and energy densities.

At $`5/2<\tau <7/2`$ the baryonic density of the gas in the mixed phase, $`\rho _g^{mix}\rho _{id}^{mix}(T)/(1+b\rho _{id}^{mix}(T))`$, becomes equal to that of the liquid at $`T>T_c`$, i.e., $`\rho _g^{mix}=1/b\rho _\mathrm{o}`$, since

$`\rho _{id}^{mix}(T)\rho _{id}(T,\mu ^{}(T))`$ $`=`$ $`\left({\displaystyle \frac{mT}{2\pi }}\right)^{3/2}[z_1\mathrm{exp}({\displaystyle \frac{W}{T}})`$ (10)
$`+`$ $`{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}k^{\frac{5}{2}\tau }\mathrm{exp}({\displaystyle \frac{\sigma k^{2/3}}{T}})]\mathrm{},`$
if surface tension vanishes $`\sigma =0`$. It is easy to prove that the entropy and energy densities for the liquid and gas phases are also equal to each other. There are discontinuities only in the derivatives of these densities over $`T`$ and $`\mu `$, i.e., $`p(T,\mu )`$ has discontinuities of its second derivatives. Therefore, the liquid-gas transition at $`T>T_c`$ for $`5/2<\tau <7/2`$ becomes the 2-nd order phase transition. According to standard definition, the point $`T=T_c`$, $`\rho =1/b`$ separating the first and second order transitions is the tricritical point. One can see that this point is now at a finite pressure.

It is interesting to note that at $`\tau >0`$ the mixed phase boundary shown in Fig.5 is not so steep function of $`T`$ as in the case $`\tau =0`$ presented in Fig.1. Therefore, the peak in the specific heat discussed above becomes less pronounced.

In conclusion, the simplified version of the SMM is solved analytically in the grand canonical ensemble. The progress is achieved by reducing the description of phase transitions to the investigation of the isobaric partition function singularities. The model clearly demonstrates a 1-st order phase transition of the liquid-gas type. The considered system has peculiar properties near the boundary of the mixed and gaseous phases. The rapid change of the thermodynamical functions with $`T`$ at fixed $`\rho `$ takes place near this boundary due to the disappearance of the infinite liquid fragment. This leads to leveling of the caloric curves shown in Fig. 3 at temperatures between 6 – 10 MeV depending on the density. As a consequence a sharp peak and a finite discontinuity are developed in the specific heat $`c_\rho (T)`$ at the boundary of the mixed and gaseous phases.

The phase diagram appears to be rather sensitive to the value of the parameter $`\tau `$ in the Fisher’s free energy term included in our treatment. New interesting possibilities for the phase diagram emerge for $`\tau >5/2`$ in comparison with the case when $`\tau <5/2`$. The case $`5/2<\tau <7/2`$ is particularly interesting because of the appearance of the tricritical point separating the 1-st and 2-nd order phase transitions.

Acknowledgments.

The authors are grateful to D. Blaschke, J. Polonyi and P.T. Reuter for stimulating discussions. We thank the Alexander von Humboldt Foundation and DAAD (Germany) for the financial support. The research described in this publication was made possible in part by Award No. UP1-2119 of the U.S. Civilian Research & Development Foundation for the Independent States of the Former Soviet Union (CRDF). K.A.B. gratefully acknowledges the warm hospitality of the Particle- and Astrophysics Group of the University of Rostock, where this talk was given.


format: text

text:
# Tensor-Scalar Torsion

## 1 Introduction

Recently, a theory of gravitation with torsion of the form111 Here, and in the remainder of this paper, we shall denote partial differentiation by comma and (anti)symmetrization over multiple indices by (square brackets) parentheses.

$$S_{\mu \nu \sigma }=\psi _{[\mu \nu ,\sigma ]},$$
(1)
where $`\psi _{\mu \nu }`$ is an antisymmetric torsion potential has been investigated . There are several reasons for our interest in this approach. The theory has predictive power owing to the fact that second-order differential equations for $`\psi _{\mu \nu }`$ are obtained when the curvature scalar is used as the Lagrangian . This theory makes definite predictions of a long range spin-spin interaction , helicity flip interactions , and radiation of torsion waves , although the effects may be smaller than can be measured with present technology. In this theory, the source for torsion is intrinsic spin, and it has been shown that the correct law for conservation of rotational angular momentum plus spin is consistently obtained.

One motivation for analysis of gravitational theories with torsion is the generic appearance of the antisymmetric tensor field in models which correspond to the low-energy limit of string theory. The argument that (1) may be related to string theory first requires the expansion of the curvature scalar $`R`$ for $`U_4`$ spacetime222 In a $`U_n`$ spacetime the metric or fundamental tensor is real-valued and symmetric, and the connection is metrically compatible. The additional constraint that the connection be symmetric defines a $`V_n`$ spacetime. The details are to be found in . Our notation is such that covariant differentiation with the Christoffel connection formed from the metric, $`\{_{\alpha \beta }^\mu \}`$, shall be indicated by a semi-colon, while the full covariant derivative, using the entire affine connection including torsion, is expressed by the operator $`_\mu `$. in terms of the $`V_4`$ scalar curvature $`{}_{}{}^{o}R`$ plus torsion terms

$$R={}_{}{}^{o}RS^{\alpha \beta \gamma }S_{\alpha \beta \gamma }.$$
(2)
This is equivalent to the low energy string theory Lagrangian, for a constant dilaton .

Another reassuring feature of this theory is that when matter sources are introduced in the usual manner, consistent physical results are obtained for phenomenological spin sources and for classical Dirac fields which possess an intrinsic spin . However, in solving the torsion field equations and the Dirac equation under the ansatz (1), a scalar field had to be introduced as a function of integration. Denoting this additional field by $`\chi `$, the solution for the torsion $`S^{\mu \nu \sigma }`$ in the presence of a source field $`\mathrm{\Psi }`$ which obeys the Dirac equation is

$$S^{\mu \nu \sigma }=\frac{2\pi i\mathrm{}G}{c^3}\overline{\mathrm{\Psi }}\gamma ^{[\mu }\gamma ^\nu \gamma ^{\sigma ]}\mathrm{\Psi }+ϵ^{\mu \nu \sigma \alpha }\chi _{,\alpha }.$$
(3)
Introducing the dual torsion $`b_\mu `$, defined as $`b_\mu =ϵ_{\mu \alpha \beta \gamma }S^{\alpha \beta \gamma }`$, this may be rewritten as

$$b^\mu =\frac{12\pi \mathrm{}G}{c^3}\overline{\mathrm{\Psi }}\gamma ^\mu \gamma _5\mathrm{\Psi }6\chi ^{,\mu }.$$
(4)
To see that the presence of the scalar field is vital to the dynamics of the general solution, consider the following argument. From the ansatz (1), it follows that $`b^\sigma {}_{;\sigma }{}^{}=0`$ must hold identically. For the $`b^\mu `$ in (4), this leads to

$$\mathrm{}\chi \chi ^{,\mu }{}_{;\mu }{}^{}=\frac{12\pi \mathrm{}G}{c^3}\left(\overline{\mathrm{\Psi }}\gamma ^\mu \gamma _5\mathrm{\Psi }\right)_{;\mu },$$
(5)
the right hand side of which, in general, is not zero. Therefore, the scalar field is dynamical and may not be consistently chosen to vanish.

In the present paper, we investigate the consequences of naturally incorporating a scalar field a priori by modifying (1). This seems reasonable on two grounds. The first is that, as shown above, a scalar field naturally arises in the solutions of this model. The second is that the low-energy limit of string theory possesses scalar modes in addition to the tensor and vector modes. Perhaps one of these modes contributes to the gravitational sector in the way that we propose. Thus, we shall assume from the start that the torsion tensor admits the decomposition

$$S_{\mu \nu \sigma }=H_{\mu \nu \sigma }+ϵ_{\mu \nu \sigma \alpha }\lambda ^{,\alpha },$$
(6)
where $`H_{\mu \nu \sigma }=\psi _{[\mu \nu ,\sigma ]}`$, and $`\lambda `$ is the scalar field.

The main result that we obtain is that no harm is done to the physically desirable law of conservation of total angular momentum (rotation plus spin). Furthermore, it is still the case in our model that spin acts as the source for the $`H_{\mu \nu \sigma }`$ part of the torsion, while the scalar field does not directly arise from spin.

## 2 Field equations

The action integral is given by

$$I=(\sqrt{g}R+k_M)d^4x,$$
(7)
where $`\sqrt{g}R`$ and $`_M`$ are the gravitational and material Lagrangian densities respectively, and $`k=8\pi G/c^4`$. The curvature scalar, $`R=g^{\alpha \beta }R_{\alpha \beta }`$, includes torsion terms. With torsion of the form (6), $`R`$ expands as

$$R={}_{}{}^{o}RH^{\alpha \beta \gamma }H_{\alpha \beta \gamma }+2ϵ^{\alpha \beta \gamma \sigma }H_{\beta \gamma \sigma }\lambda _{,\alpha }+6\lambda ^{,\alpha }\lambda _{,\alpha }.$$
(8)
Variations of $`I`$ with respect to $`g_{\mu \nu }`$, $`\psi _{\mu \nu }`$, and $`\lambda `$, are assumed to vanish, yielding the classical field equations. With the following notational conventions for the sources of the metric, antisymmetric potential, and scalar fields:

$$\sqrt{g}T^{(\mu \nu )}\delta g_{\mu \nu }d^4x\frac{\delta _M}{\delta g^{\mu \nu }}\delta g_{\mu \nu }d^4x,$$
(9)
$$\sqrt{g}T^{[\mu \nu ]}\delta \psi _{\mu \nu }d^4x\frac{\delta _M}{\delta \psi ^{\mu \nu }}\delta \psi _{\mu \nu }d^4x,$$
(10)
$$\sqrt{g}U\delta \lambda d^4x\frac{\delta _M}{\delta \lambda }\delta \lambda d^4x,$$
(11)
the field equations are as follows:

$${}_{}{}^{o}G_{}^{\mu \nu }3H^\mu {}_{\alpha \beta }{}^{}H_{}^{\nu \alpha \beta }+\frac{1}{2}g^{\mu \nu }H^{\alpha \beta \gamma }H_{\alpha \beta \gamma }+6\lambda ^{,\mu }\lambda ^{,\nu }3g^{\mu \nu }\lambda ^{,\alpha }\lambda ^{,\alpha }=kT^{(\mu \nu )},$$
(12)
$$H^{\mu \nu \sigma }{}_{;\sigma }{}^{}=\frac{k}{2}T^{[\mu \nu ]},$$
(13)
$$\mathrm{}\lambda =\frac{U}{12}.$$
(14)
Note that the use of the symbol $`T^{[\mu \nu ]}`$ to denote the source of the antisymmetric field is to allow the formal introduction of $`T^{\mu \nu }=T^{(\mu \nu )}+T^{[\mu \nu ]}`$ as the source for the combined metric and torsion potential $`g^{(\mu \nu )}+\psi ^{[\mu \nu ]}`$. $`U`$ describes the phenomenological source for the scalar field and $`{}_{}{}^{o}G`$ is the $`V_4`$ Einstein tensor.

With these field equations we shall develop the equations of motion for a test body with intrinsic spin moving through a region of spacetime containing both gravitational and torsion fields. Furthermore, we verify the consistency of the conservation law for rotational plus intrinsic spin angular momentum in the case in which the external torques due to gravity and torsion vanish.

## 3 Equations of motion for a small test body

In reference  Hammond adapted the method established by Papapetrou  to develop the equations of motion for a small body in a gravitational and torsion field (no scalar field). This method is used here to write the equations of motion for a small body in a combined gravitational and torsion field in which the scalar field is present.

Papapetrou solved for the motion of a spinning test particle in a background gravitational field through the effective use of an iterative multipole-like expansion. There are systematic and practical advantages to this approach . Spinless and structureless extended particles are zeroth-order poles, while spinning particles possess at least a pole–dipole structure. For small bodies it is possible to neglect higher-order pole terms to a reasonable approximation.

The $`U_4`$ Bianchi Identity  holds when torsion is present.

$$_\nu G^{\mu \nu }=2S^{\mu \alpha \beta }R_{\beta \alpha }S_{\alpha \beta \gamma }R^{\mu \gamma \beta \alpha }.$$
(15)
Substitution of the field equations and use of $`U_4`$ geometrical identities leads to

$$T^{\mu \nu }{}_{;\nu }{}^{}=\frac{3}{2}H^\mu {}_{\alpha \beta }{}^{}T_{}^{\alpha \beta }+\frac{1}{6k}\lambda ^{,\mu }U,$$
(16)
where $`T^{\mu \nu }=T^{(\mu \nu )}+T^{[\mu \nu ]}`$, as was discussed above.

Let $`\tau ^{\mu \nu }=\sqrt{g}T^{(\mu \nu )}`$, $`j^{\alpha \beta }=\frac{k}{K}T^{[\alpha \beta ]}`$, and $`\stackrel{~}{j}^{\alpha \beta }=\sqrt{g}j^{\alpha \beta }`$ so that a tilde denotes density, where $`K`$ is an arbitrary coupling constant for the spin field introduced by Hammond in . It then follows that

$$\tau ^{\mu \nu }{}_{,\nu }{}^{}=\sqrt{g}T^{(\mu \nu )}{}_{;\nu }{}^{}\{_{\alpha \beta }^\mu \}\tau ^{\alpha \beta },$$
(17)
and therefore

$$\tau ^{\mu \nu }{}_{,\nu }{}^{}=\frac{3K}{2k}H^\mu {}_{\alpha \beta }{}^{}\stackrel{~}{j}_{}^{\alpha \beta }+\frac{1}{6k}\lambda ^{,\mu }\stackrel{~}{U}\{_{\alpha \beta }^\mu \}\tau ^{\alpha \beta }.$$
(18)
Consider the integral of $`\tau ^{\mu \nu }_{,\nu }`$ over a spatial hypersurface characterized by fixed time. For a bounded source, the integral has support only over the finite volume of the body, and hence the surface terms which arise upon invoking the divergence theorem vanish (viz., $`\tau ^{\mu i}{}_{,i}{}^{}dV=0`$, where Latin indices run over spatial coordinates). Then,

$$\tau ^{\mu 0}{}_{,0}{}^{}dV=\{\frac{3K}{2k}H^\mu {}_{\alpha \beta }{}^{}\stackrel{~}{j}_{}^{\alpha \beta }+\frac{1}{6k}\lambda ^{,\mu }\stackrel{~}{U}\{_{\alpha \beta }^\mu \}\tau ^{\alpha \beta }\}dV.$$
(19)
Let $`y^\alpha `$ denote the position of the center of mass of the small test body, while $`\delta x^\alpha `$ are the center of mass coordinates of points in the body. Thus $`x^\alpha =y^\alpha +\delta x^\alpha `$ describes the position of any point within the body. Expanding $`\{_{\alpha \beta }^\mu \}`$, $`H^\mu _{\alpha \beta }`$, and $`\lambda ^{,\mu }`$ in Taylor series about $`y^\alpha `$ to first order in the small quantities $`\delta x^\alpha `$, and substituting these into (19) results in

$`{\displaystyle \frac{d}{dt}}{\displaystyle \tau ^{\mu 0}𝑑V}`$ $`=`$ $`{\displaystyle }\{{\displaystyle \frac{3K}{2k}}(H^\mu {}_{\alpha \beta }{}^{}+H^\mu {}_{\alpha \beta ,\sigma }{}^{}\delta x^\sigma )\stackrel{~}{j}^{\alpha \beta }`$
$`+{\displaystyle \frac{1}{6k}}(\lambda ^{,\mu }+\lambda ^{,\mu }{}_{,\sigma }{}^{}\delta x^\sigma )\stackrel{~}{U}`$
$`(\{_{\alpha \beta }^\mu \}+\{_{\alpha \beta }^\mu \}_{,\sigma }\delta x^\sigma )\tau ^{\alpha \beta }\}dV.`$
The time derivative has been pulled out of the integral on the left hand side of (3), and on the right hand side it is understood that the fields and their derivatives are evaluated at the position of the center of mass.

Equation (3) is the equation of motion for the body. To see this more clearly, we note that $`v^\alpha =dy^\alpha /d\tau `$ is the center of mass velocity of the body ($`\tau `$ is the invariant proper time as measured along the worldline of the center of mass), and define the following four quantities:

$$M^{\mu \nu }\frac{v^0}{c}\tau ^{\mu \nu }𝑑V,$$
(21)
$$M^{\alpha \mu \nu }\frac{v^0}{c}\delta x^\alpha \tau ^{\mu \nu }𝑑V,$$
(22)
$$J^{\mu \nu }\frac{1}{c}(\delta x^\mu \tau ^{\nu 0}\delta x^\nu \tau ^{\mu 0})𝑑V,$$
(23)
$$m^{\alpha \mu \nu }\frac{v^0}{c}\delta x^\alpha \stackrel{~}{j}^{\mu \nu }𝑑V.$$
(24)
In the language used by Papapetrou, $`M^{\mu \nu }`$ consists of pole terms, while $`M^{\alpha \mu \nu }`$ and $`J^{\mu \nu }`$ are dipole contributions. $`m^{\alpha \mu \nu }`$ is the lowest-order torsion term in this approximation scheme, and is also a dipole term. Reexpressing (3) in terms of these quantities produces

$`{\displaystyle \frac{d}{d\tau }}{\displaystyle \frac{M^{\mu 0}}{v^0}}`$ $`=`$ $`{\displaystyle \frac{3K}{2k}}{\displaystyle \frac{v^0}{c}}H^\mu {\displaystyle {}_{\alpha \beta }{}^{}}\stackrel{~}{j}^{\alpha \beta }𝑑V+{\displaystyle \frac{3K}{2k}}H^\mu {}_{\alpha \beta ,\gamma }{}^{}m_{}^{\gamma \alpha \beta }`$
$`\{_{\alpha \beta }^\mu \}M^{\alpha \beta }+\{_{\alpha \beta }^\mu \}_{,\gamma }M^{\gamma \alpha \beta }`$
$`+{\displaystyle \frac{1}{6}}{\displaystyle \frac{v^0}{ck}}{\displaystyle \lambda ^{,\mu }\stackrel{~}{U}𝑑V}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{v^0}{ck}}{\displaystyle \delta x^\gamma \lambda ^{,\mu }{}_{,\gamma }{}^{}\stackrel{~}{U}𝑑V},`$
which, if $`H^\mu {}_{\alpha \beta }{}^{}=0`$ and $`\lambda =0`$, reduces to the familiar set of equations for a rotating body in a gravitational field. Equation (3) is interpreted as follows, the term on the left hand side is the net force acting on the test body, and the terms on the right hand side are the forces due to the spin field and its gradient, the gravitational field and its gradient, and the scalar field and its gradient.

The equations for angular momentum and spin are obtained by integrating the divergence of the first moment of the symmetric stress-energy-momentum density,

$$(x^\sigma \tau ^{\mu \nu }),_\nu =\tau ^{\mu \sigma }+x^\sigma \frac{3K}{2k}H^\mu {}_{\alpha \beta }{}^{}\stackrel{~}{j}_{}^{\alpha \beta }x^\sigma \{_{\alpha \beta }^\mu \}\tau ^{\alpha \beta }+\frac{1}{12k}x^\sigma \lambda ^{,\mu }\stackrel{~}{U},$$
(26)
over all space. The boundedness of the source yields $`(x^\mu \tau ^{\nu i})_{,i}𝑑V=0`$, simplifying the integral of the left hand side of (26).

$`{\displaystyle \frac{d}{dt}}{\displaystyle (x^\sigma \tau ^{\mu 0})𝑑V}`$ $`=`$ $`{\displaystyle \tau ^{\mu \sigma }𝑑V}+{\displaystyle x^\sigma \frac{3K}{2k}H^\mu {}_{\alpha \beta }{}^{}\stackrel{~}{j}_{}^{\alpha \beta }𝑑V}`$
$`{\displaystyle }x^\sigma \{_{\alpha \beta }^\mu \}\tau ^{\alpha \beta }dV+{\displaystyle }{\displaystyle \frac{1}{12k}}x^\sigma \lambda ^{,\mu }\stackrel{~}{U}dV.`$
Once again, expanding $`\{_{\alpha \beta }^\mu \}`$, $`H^\mu _{\alpha \beta }`$, and $`\lambda ^{,\mu }`$ about $`y^\alpha `$, and using the quantities defined above, we see that

$`M^{\mu \sigma }`$ $`=`$ $`{\displaystyle \frac{v^\sigma }{v^0}}M^{\mu 0}{\displaystyle \frac{d}{d\tau }}{\displaystyle \frac{M^{\sigma \mu 0}}{v^0}}{\displaystyle \frac{3K}{2k}}H^\mu {}_{\alpha \beta }{}^{}m_{}^{\sigma \alpha \beta }`$
$`\{_{\alpha \beta }^\mu \}M^{\sigma \alpha \beta }{\displaystyle \frac{v^0}{12kc}}\lambda ^{,\mu }{\displaystyle }\delta x^\sigma \stackrel{~}{U}dV.`$
Following upon its construction, $`M^{\mu \sigma }`$ must be symmetric in its indices, which requires, after some simplification, that

$`{\displaystyle \frac{d}{d\tau }}J^{\sigma \mu }+v^\sigma {\displaystyle \frac{M^{\mu 0}}{v^0}}v^\mu {\displaystyle \frac{M^{\sigma 0}}{v^0}}`$ $`=`$ $`{\displaystyle \frac{3K}{k}}H^{[\mu }{}_{\alpha \beta }{}^{}m_{}^{\sigma ]\alpha \beta }+2\{_{\alpha \beta }^{[\mu }\}M^{\sigma ]\alpha \beta }`$
$`+{\displaystyle \frac{v^0}{kc}}\lambda ^{,[\mu }{\displaystyle \delta x^{\sigma ]}\stackrel{~}{U}𝑑V}.`$
The right hand side of (3) is the net applied torque on the body due to the spin field, the gravitational field, and the scalar field. The terms on the left side represent the proper time derivatives of the rotational angular momentum plus spin, and the orbital angular momentum.

As defined in (23), $`J^{\sigma \mu }`$ is the total angular momentum. In Hammond’s original model and again in the present case it is composed of $`L^{\sigma \mu }`$, the rotational angular momentum, plus $`S^{\sigma \mu }`$, arising from intrinsic spin which is not associated with any physical motion of the body:

$$J^{\sigma \mu }L^{\sigma \mu }+S^{\sigma \mu }.$$
(30)
When the applied net torque from all sources vanishes, (3) may be integrated along the worldline of the test body to yield an expression of the law of conservation of angular momentum

$$L^{\sigma \mu }+S^{\sigma \mu }+y^\sigma p^\mu p^\sigma y^\mu =\mathrm{constant}.$$
(31)
## 4 Conclusion

The existence of rather general arguments in favor of the addition of a scalar field to the torsion sector of Hammond’s theory was in no way a guarantee that such an augmented model would remain self-consistent. In this paper we have shown that the inclusion of the scalar field in the specified manner with a classical phenomenological source may be achieved without a breakdown in the conservation laws governing the motion of classical test bodies. This result suggests that the model deserves further investigation.

The only serious deficiency that we have noted is that the source of the scalar field is nearly unconstrained classically – unlike the case of the $`\psi _{[\mu \nu ]}`$ torsion potential in Hammond’s original model. Now that the augmented model has passed its first hurdle, we plan to incorporate a classical Dirac field as the source of stress-energy and spin. It is hoped that this analysis may shed light upon the physical origin and properties of the additional scalar field.

## 5 Acknowledgements

This research was funded under a Department Of Energy Traineeship, NASA grant NAG81007, NSF EPSCoR grant OSR-9452892, and a North Dakota State University Grant-In-Aid.


text:
# An Analysis of the Quasicontinuum Method

### Abstract

The aim of this paper is to present a streamlined and fully three-dimensional version of the quasicontinuum (QC) theory of Tadmor et al.  and to analyze its accuracy and convergence characteristics. Specifically, we assess the effect of the summation rules on accuracy; we determine the rate of convergence of the method in the presence of strong singularities, such as point loads; and we assess the effect of the refinement tolerance, which controls the rate at which new nodes are inserted in the model, on the development of dislocation microstructures.

Keywords: A. Quasicontinuum; B. Error analysis; C. Nanoindentation; D. Atomistic models

## 1 Introduction

The aim of this paper is to present a streamlined and fully three-dimensional version of the quasicontinuum (QC) theory of Tadmor et al.  and to analyze its accuracy and convergence characteristics. The theory of the quasicontinuum furnishes a computational scheme for seamlessly bridging the atomistic and continuum realms. The chief objective of the theory is to systematically coarsen an atomistic description by–and only by–the judicious introduction of kinematic constraints. These kinematic constraints are selected and designed so as to preserve full atomistic resolution where required, e. g., in the vicinity of lattice defects, and to treat collectively large numbers of atoms in regions where the deformation field varies slowly on the scale of the lattice. Thus, in its purest form all input into the theory concerning the behavior of the material is atomistic, and all approximations are strictly kinematic in nature.

Different variants of the theory have been developed and documented over a series of publications , where numerous examples of application have also been presented. Here we endeavor to emphasize the essential building blocks of the static theory, to wit: i) the constrained minimization of the atomistic energy of the solid; ii) the use of summation rules in order to compute the effective equilibrium equations; iii) and the use of adaption criteria in order to tailor the computational mesh to the structure of the deformation field. In particular, we develop a new class of summation rules based on sampling lattice functions over clusters of atoms. In Section 4, we present a detailed numerical analysis that probes various aspects of the accuracy and performance of the method. Specifically, we assess the effect of the cluster size and lumping procedure on accuracy; we determine the rate of convergence of the method in the presence of strong singularities, such as point loads; and we assess the effect of the refinement tolerance, which controls the rate at which new nodes are inserted in the model, on the development of dislocation microstructures.

## 2 The quasicontinuum method

We consider a reference configuration of a crystal such that its $`N`$ atoms occupy a subset of a simple $`d`$-dimensional Bravais lattice spanned by basis vectors $`\left\{𝒂_i,i=1,\mathrm{},d\right\}`$. The coordinates of the atoms in this reference configuration are:

$$𝑿\left(𝒍\right)=\underset{i=1}{\overset{d}{}}l^i𝒂_i,𝒍Z^d.$$
(1)
Here, $`𝒍`$ are the lattice coordinates which designate individual atoms, $`Z`$ is the set of integer numbers and $`d`$ is the dimension of space. From the definition, it follows that the set $``$ is the collection of lattice sites occupied by atoms. The coordinates of the atoms in a deformed configuration of the crystal are denoted $`\left\{𝒒\left(𝒍\right),𝒍\right\}`$. For convenience, we shall collect all atomic coordinates in an array $`𝒒`$ and regard such array as an element of the linear space $`XR^{Nd}`$, the ‘configuration’ space of the crystal.

The energy of the crystal is assumed to be expressible as a function $`E\left(𝒒\right)`$, e. g., through the use of empirical interatomic potentials. In addition, the crystals may be subjected to applied loads. We assume these loads to be conservative and to derive from an external potential $`\mathrm{\Phi }^{\mathrm{ext}}\left(𝒒\right)`$. The total potential energy is, therefore:

$$\mathrm{\Phi }\left(𝒒\right)=E\left(𝒒\right)+\mathrm{\Phi }^{\mathrm{ext}}\left(𝒒\right)$$
(2)
The crystal may also be subject to displacement boundary conditions over part of its boundary, e. g., as a result of the application of a rigid indentor. We wish to determine the stable equilibrium configurations of the crystal. The problem is, therefore, to determine all the *local* minima of $`\mathrm{\Phi }\left(𝒒\right)`$ consistent with the displacement, or essential, boundary conditions. We shall, somewhat equivocally, enunciate this problem as:

$$\underset{𝒒X}{\mathrm{min}}\mathrm{\Phi }\left(𝒒\right)$$
(3)
We emphasize, however, that in general the aim is not simply to determine the absolute minimizer of $`\mathrm{\Phi }\left(𝒒\right)`$, but rather the much richer, and physically more relevant, set of metastable configurations of the crystal.

The distinction arises from the fact that the energy function $`E\left(𝒒\right)`$ must be invariant under the class of affine deformations which map the crystal lattice onto itself. This implies, for instance, that $`E\left(𝒒\right)`$ must be periodic under crystallographic slip with period equal to the Burgers vector. Hence, the energy function $`E\left(𝒒\right)`$ is strongly nonconvex and the potential energy $`\mathrm{\Phi }\left(𝒒\right)`$ in general has vast numbers of local minimizers, or metastable configurations. Techniques for systematically exploring the full set of metastable configurations, and indeed the entire configuration space of the crystal, include statistical mechanics and path-integral methods. An alternative approach is to load the crystal incrementally and proceed by continuation. Thus, at every step in the solution procedure the loads, or the prescribed displacements, are incremented and the system is allowed to relax to a nearby stable configuration. The relaxation process may be governed by inertia, viscosity, kinetics, or some other rate-limiting mechanism; or may be numerical in nature and owe to the use of iterative solvers such as conjugate gradients. This latter approach is adopted in all the numerical tests presented in this paper.

### 2.1 Interpolation

The essence of the theory of Tadmor et al.  is to replace (3) by a constrained minimization of $`\mathrm{\Phi }\left(𝒒\right)`$ over a suitably chosen subspace $`X_h`$ of $`X`$. In order to define $`X_h`$, we begin by selecting a reduced set $`_h`$ of $`N_h<N`$ ‘representative atoms’ (Fig. 1). The selection of the representative atoms is based on the local variation of the fields and will be addressed subsequently. In addition, we introduce a triangulation $`𝒯_h`$ of $`_h`$. It bears emphasis that the triangulation $`𝒯_h`$ may be unstructured. In particular, $`_h`$ need not define a Bravais lattice. The positions of the remaining atoms are determined by piecewise linear interpolation of the representative atom coordinates. We shall regard the resulting coordinates $`\left\{𝒒_h\left(𝒍\right),𝒍\right\}`$ as belonging to a linear space $`X_h`$ of dimension $`N_hd`$.

Let $`\phi _h\left(𝑿|𝒍_h\right)`$, $`𝒍_h_h`$, be a collection of shape functions for $`𝒯_h`$. Thus, $`\phi _h\left(𝑿|𝒍_h\right)`$ is continuous and piecewise linear, its domain is restricted to the simplices $`K𝒯_h`$ incident to $`𝑿\left(𝒍_h\right)`$, and it vanishes at all nodes of the triangulation except at $`𝑿\left(𝒍_h\right)`$, where it takes the value $`1`$, i. e.,

$$\phi _h\left(𝑿\left(𝒍_h^{}\right)|𝒍_h\right)=\delta \left(𝒍_h^{}|𝒍_h\right).$$
(4)
By construction,

$$𝒒_h\left(𝒍\right)=\underset{𝒍_h_h}{}\phi _h\left(𝒍|𝒍_h\right)𝒒_h\left(𝒍_h\right),$$
(5)
where we write

$$\phi _h\left(𝒍|𝒍_h\right)=\phi _h\left(𝑿\left(𝒍\right)|𝒍_h\right).$$
(6)
Evidently, $`\left\{\phi _h\left(𝒍|𝒍_h\right),𝒍_h_h\right\}`$ constitutes a basis for $`X_h`$ and the fields $`𝒒_h\left(𝒍\right)`$ are entirely determined by their values $`𝒒_h\left(𝒍_h\right)`$ at the representative atoms. In addition, the basis lattice functions are required to satisfy the identity:

$$\underset{𝒍_h_h}{}\phi _h\left(𝒍|𝒍_h\right)=1$$
(7)
i. e., the basis lattice functions must define a partition of unity over $``$. This requirement ensures that constant fields are interpolated exactly by the basis lattice functions.

### 2.2 Reduced equations

The reduced counterpart of problem (3) is now

$$\underset{𝒒_hX_h}{\mathrm{min}}\mathrm{\Phi }\left(𝒒_h\right).$$
(8)
The minimizers of the reduced problem satisfy the reduced equations of equilibrium

$$𝒇_h\left(𝒍_h\right)=\underset{𝒍}{}𝒇\left(𝒍|𝒒_h\right)\phi _h\left(𝒍|𝒍_h\right)=\mathrm{𝟎}.$$
(9)
Here,

$$𝒇\left(𝒒\right)=\mathrm{\Phi },_𝒒\left(𝒒\right)$$
(10)
are the forces corresponding to $`𝒒`$ and $`𝒇\left(𝒍|𝒒\right)`$ is the value of $`𝒇\left(𝒒\right)`$ at site $`𝒍`$. Thus, the reduced problem entails the solution of the $`N_hd`$ equations (9) in the $`N_hd`$ unknowns $`𝒒\left(𝒍_h\right)`$, $`𝒍_h_h`$.

### 2.3 Node-based summation rules

The practicality of the method further hinges on the possibility of avoiding the calculation of the full atomistic force array $`𝒇`$ and carrying out full lattice sums, as seemingly required in (9). As noted by Tadmor et al. , this may be accomplished by the introduction of summation rules similar to the conventional quadrature rules of numerical integration. The problem is thus to approximate sums of the general form

$$S=\underset{𝒍}{}g\left(𝒍\right),$$
(11)
where $`g\left(𝒍\right)`$ is a lattice function. By analogy to numerical quadrature rules, we begin by exploring summation rules of the form

$$S\underset{𝒍𝒮_h}{}n_h\left(𝒍\right)g\left(𝒍\right)S_h$$
(12)
for some suitably chosen collection of summation points $`𝒮_h`$ and weights $`n_h\left(𝒍\right)`$, not necessarily integer. Loosely speaking, the weights $`n_h\left(𝒍\right)`$ may be regarded as the number of atoms represented by the site $`𝒍`$. Proceeding as in the development of numerical quadrature rules, the weights $`n_h\left(𝒍\right)`$ may be determined by the requirement that the summation rule (12) be exact for a restricted class of lattice functions .

The lowest-order summation rule is obtained by requiring that all lattice functions in $`X_h`$, i. e., all piecewise linear functions supported by the triangulation $`𝒯_h`$, be summed exactly. This is tantamount to setting $`𝒮_h=_h`$ and requiring that the summation rule (12) be exact for all shape functions $`\phi _h\left(𝒍|𝒍_h\right)`$, $`𝒍_h_h`$. This requirement gives, explicitly,

$$n_h\left(𝒍_h\right)=\underset{𝒍}{}\phi _h\left(𝒍|𝒍_h\right),𝒍_h_h.$$
(13)
In fine regions of the triangulation, the sums on the right-hand side of (13) may be computed explicitly. By contrast, in coarse regions of the mesh approaching the continuum limit this explicit calculation becomes impractical. However, in such regions the lattice sum (12) ostensibly reduces to an integral and the corresponding weights are those of conventional Lobatto quadrature . In this limit, each simplex $`K𝒯_h`$ simply contributes $`\left(N/V\right)\left|K\right|/\left(d+1\right)`$ atoms to each of its $`d+1`$ nodes, where $`N/V`$ is the atom density in the undeformed configuration of the crystal.

As a simple illustrative example, consider a monoatomic chain discretized into ‘elements’ each containing $`L`$ bonds. For an internal node, the summation weight follows as

$$n_h=1+2\underset{l=1}{\overset{L1}{}}\left(1\frac{l}{L}\right)=L,$$
(14)
which coincides with the continuum limit. By contrast, for an end node the summation weight is

$$n_h=1+\underset{l=1}{\overset{L1}{}}\left(1\frac{l}{L}\right)=\frac{L+1}{2}.$$
(15)
We note that $`n_h`$ tends to the continuum Lobatto limit of $`L/2`$ as $`L\mathrm{}`$ and reduces to $`n_h=1`$ in the full atomistic limit, as required.

### 2.4 Reduced equations based on summation rules

The application of a summation rule to the evaluation of the equilibrium equations (9) leads to the system of equations

$$𝒇_h\left(𝒍_h\right)\underset{𝒍𝒮_h}{}n_h\left(𝒍\right)𝒇\left(𝒍|𝒍_h\right)\phi _h\left(𝒍|𝒍_h\right)=\mathrm{𝟎}.$$
(16)
As desired, the calculation of the effective forces in (16) is of complexity $`𝒪\left(N_h\right)`$ provided that the number of sampling sites in $`𝒮_h`$ is of order $`N_h`$ and the atomic interactions are short-ranged.

The formulation of the method takes a somewhat unexpected turn at this point in that the node-based summation rule discussed in the foregoing suffers from a rank-deficiency problem. Thus, in the linear case the node-based summation rule leads to rank-deficient systems of equations. This is in analogy to the rank-deficiency of finite-element stiffness matrices computed by Lobatto quadrature, and manifests itself as the existence of so-called zero-energy modes which pollute the displacement field or render the system of equations singular altogether.

By way of illustration, we consider a Lennard-Jones crystal in the form of a cube containing $`32\times 32\times 32`$ fcc unit cells. In this example, eight representative atoms are introduced at the corners of the cube, and the cube is triangulated as shown in Fig. 2a. The stability of the crystal in its reference configuration may be assessed by examining the spectral properties of the reduced stiffness matrix. A direct calculation of its eigenvalues reveals that seven of them are identically equal to zero, which in turn denotes the presence of one zero-energy deformation mode, Fig. 2b. However, the zero-energy deformation mode is absent when the lattice sums are carried out in full, i. e., when the test is repeated with the equilibrium equations (16) replaced by (9).

The numerical test just described shows that the node-based summation rules are indeed rank-deficient, but it also suggests that the reduced system is stable provided that a sufficient number of sampling points is included in the summation rules. A new set of summation rules which offers this added flexibility is presented next.

### 2.5 Cluster summation rules

A class of summation rules which generalizes the node-based rule may be obtained by sampling the lattice function over neighborhoods of the representative atoms. We shall refer to these neighborhoods as clusters, and the resulting summation rules as cluster summation rules. Each cluster may be regarded as a representative crystallite where the state and the behavior of the crystal are sampled.

More specifically, let $`𝒞\left(𝒍_h\right)=\{𝒍:\left|𝑿\left(𝒍\right)𝑿\left(𝒍_h\right)\right|r\left(𝒍_h\right)\}`$ be the cluster of lattice sites located within a sphere of radius $`r\left(𝒍_h\right)`$ centered on the representative atom $`𝒍_h`$ (Fig. 3). Assume for now that the clusters thus defined do not overlap. The cluster summation rule induced by the set of clusters $`𝒞\left(𝒍_h\right),𝒍_h_h`$ is

$$S_h=\underset{𝒍_h_h}{}n_h\left(𝒍_h\right)S\left(𝒍_h\right),$$
(17)
where $`S\left(𝒍_h\right)`$ denotes the sum over all atoms in the cluster $`𝒞\left(𝒍_h\right)`$, i. e.,

$$S\left(𝒍_h\right)=\underset{𝒍𝒞\left(𝒍_h\right)}{}g\left(𝒍\right).$$
(18)
As before, the cluster weights $`n_h\left(𝒍_h\right),𝒍_h_h`$ are computed by requiring that the cluster summation rule (17) be exact for all basis functions.

It should be carefully noted that, as the triangulation size approaches the atomic length scale, the spherical clusters defined previously in general overlap, Fig. 4. In this case, we truncate the clusters as follows. A site $`𝒍`$ belongs to cluster $`𝒞\left(𝒍_h\right)`$ if its distance to $`𝑿\left(𝒍_h\right)`$ is less than a prescribed value $`r\left(𝒍_h\right)`$ and less than its distance to any other representative atom. Thus, the cluster $`𝒞\left(𝒍_h\right)`$ is the intersection of the Voronoi cell containing $`𝑿\left(𝒍_h\right)`$ and sphere of radius $`r\left(𝒍_h\right)`$ centered at $`𝑿\left(𝒍_h\right)`$. Ambiguous cases corresponding to sites which are equidistant from two representative atoms are resolved randomly. In the limit of full atomistic resolution, each cluster contains exactly one atom and the corresponding weight is $`1`$. Thus, the cluster summation rules reduce seamlessly to the exact lattice sum when the fully atomistic limit is attained.

#### 2.5.1 Summation-rule error analysis for monatomic chains

The accuracy of cluster summation rules can be analyzed simply in the case of a one-dimensional monatomic chain. In this case, $`Z`$ is the set of all sites in the chain. Without loss of generality, we take the coordinate of a site $`l`$ in the undeformed configuration of the chain to be $`X\left(l\right)=l`$. For simplicity, we consider a uniform triangulation of $``$ of size $`h`$. Thus, the coordinates of the representative atoms are $`X\left(l_h\right)=l_hh`$, $`l_hZ`$, where $`h`$ is the simplex size. We consider clusters of uniform radius $`r<h/2`$, Fig. 5.

The summation weights follow from the requirement that the basis functions $`\phi _h\left(l|l_h\right)`$ be summed exactly, with the result,

$$n_h\left(l_h\right)=\frac{h}{1+2r},l_hZ.$$
(19)
Thus, for the infinite monatomic chain, $`n_h\left(l_h\right)`$ is simply the ratio of the number of atoms in a simplex to the number of atoms in a cluster.

In order to obtain the summation error to leading order, we consider the sum of a quadratic polynomial supported on two adjacent simplices in the reduced chain, namely,

$$g\left(l\right)=\{\begin{array}{cc}l^2\hfill & l[h,h]\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(20)
The exact sum of $`g\left(l\right)`$ over all sites of the chain is

$$S=\underset{l}{}g\left(l\right)=\underset{l=h}{\overset{l=h}{}}l^2=\frac{h\left(1+h\right)\left(1+2h\right)}{3}.$$
(21)
The calculation of the approximate sum $`S_h`$ requires three local cluster sums $`S_h\left(h\right)`$, $`S_h\left(0\right)`$ and $`S_h\left(h\right)`$

$`S_h\left(h\right)`$ $`=`$ $`{\displaystyle \underset{l=h}{\overset{l=h+r}{}}}l^2={\displaystyle \frac{\left(1+r\right)\left(6h^2+r6hr+2r^2\right)}{6}},`$ (22)
$`S_h\left(0\right)`$ $`=`$ $`{\displaystyle \underset{l=r}{\overset{l=r}{}}}l^2={\displaystyle \frac{r\left(1+r\right)\left(1+2r\right)}{3}},`$ (23)
$`S_h\left(h\right)`$ $`=`$ $`{\displaystyle \underset{l=hr}{\overset{l=h}{}}}l^2={\displaystyle \frac{\left(1+r\right)\left(6h^2+r6hr+2r^2\right)}{6}}.`$ (24)
The use of the cluster summation rules leads to the approximate lattice sum

$$S_h=n_h\left[S_h\left(h\right)+S_h\left(0\right)+S_h\left(h\right)\right]=\frac{2h\left(1+r\right)}{3\left(1+2r\right)}\left[2r^2+r\left(13h\right)+3h^2\right].$$
(25)
The summation error is $`\left|S_hS\right|`$. An upper bound on the error is obtained by setting $`r=0`$. Then, for $`h1`$ we have $`\left|S_hS\right|_{r=0}\left(4/3\right)h^3`$, which shows that the summation rule is of third order.

The effect of the cluster size on the accuracy of the summation rule may be understood as follows. For simplicity, we consider the case of $`h1`$, whereupon (21) simplifies to

$$S\frac{2}{3}h^3$$
(26)
and (25) to

$$S_h=\frac{2\left(1+r\right)}{3\left(1+2r\right)}\left(2\xi ^23\xi +3\right)h^3$$
(27)
where we write $`\xi =r/h`$. In the regime of $`\xi 1`$, (27) simplifies further and the corresponding relative error behaves as

$$\frac{\left|S_hS\right|}{\left|S\right|}\frac{2+r}{1+2r}$$
(28)
which is a decreasing function of $`r`$. Thus, as expected, increasing $`r`$ from 0 increases the accuracy of the summation rule. In the regime of $`r1`$ and $`rh`$, the relative error behaves as

$$\frac{\left|S_hS\right|}{\left|S\right|}\frac{1}{2}\left(2\xi 1\right)\left(\xi 1\right)$$
(29)
which vanishes for $`\xi =1/2`$, or $`r=h/2`$. Thus, the exact sum is obtained when the clusters encompass the entire lattice, as required.

### 2.6 Reduced equations based on cluster summation rules

The application of the cluster based summation rule to the reduced equilibrium equations (9) yields

$$𝒇_h\left(𝒍_h\right)\underset{𝒍_h^{}𝒮_h}{}n_h\left(𝒍_h^{}\right)\left[\underset{𝒍𝒞\left(𝒍_h^{}\right)}{}𝒇\left(𝒍\right)\phi _h\left(𝒍|𝒍_h\right)\right]=\mathrm{𝟎}.$$
(30)
The computational complexity of (30) is $`𝒪\left(N_hN_r\right)`$, where $`N_r`$ is the number of lattice sites in a cluster of radius $`r`$. For $`N_r>1`$, the computational effort is in excess of that corresponding to the node-based summation rules. Thus, the optimal value of the cluster size $`r\left(𝒍_h\right)`$ is subject two opposing requirements. On one hand, the cluster size has to contain a sufficient number of atoms to ensure the appropriate accuracy and stability of the summation rule. On the other, however, such number should remain small for the method to retain its computational efficiency. The effect of the cluster size on the accuracy is investigated subsequently by way of numerical testing.

### 2.7 Adaptive selection of representative atoms

The third key component of the quasicontinuum method is the use of mesh adaption in order to tailor the computational mesh to the structure of the deformation field. Ideally, the adaptivity itself should be driven by the energetics of the system, i. e., the mesh should return the least possible potential energy for a fixed number of representative atoms. However, in order to optimize the mesh in this manner it is necessary to know the relation between the energy, or suitable energy bounds thereof, and the mesh size. Unfortunately, unlike the continuum case, an approximation theory for discrete systems such as considered here appears to be entirely lacking at present, and such energy bounds are presently unavailable.

Under these conditions, it becomes necessary to resort to empirical adaption indicators. In calculations we adopt as adaption indicator $`ϵ\left(K\right)`$ for simplex $`K`$ the variation of the displacement field over $`K`$ modulo rotations. This variation should be well-defined since the displacement field of the crystal may reasonably be expected to be of bounded variation. Specifically, we measure the variation of the displacement field as

$$ϵ\left(K\right)\sqrt{\left|II_E\left(K\right)\right|}h\left(K\right)$$
(31)
where $`II_E\left(K\right)`$ denotes the second invariant of the Lagrangian strain tensor in simplex $`K`$, and $`h\left(K\right)`$ is the size of $`K`$ (cf ). It follows from its definition that $`ϵ\left(K\right)`$ is invariant under rotations. The element $`K`$ is deemed acceptable if

$$\frac{ϵ\left(K\right)}{b}<TOL$$
(32)
for some prescribed tolerance $`TOL<1`$, and is targeted for refinement otherwise. In (32) $`b`$ denotes the magnitude of the smallest Burgers vector of the crystal.

In our implementation of the method, the refinement of the simplices which violate (32) is accomplished by the application of subdivision rules. In particular, we choose to bisect the longest edge of the simplices which are targeted for refinement. The new representative atoms are inserted in the lattice site nearest to the geometrical midpoint of the longest edge. In order to preserve the quality of the computational mesh, we apply standard mesh-improvement operations such as edge-face or octahedral swapping  and smoothing . Each new mesh is equilibrated and the remeshing criterion (32) is re-evaluated for the new solution. The process is repeated until the mesh remains unchanged.

The significance of (32) becomes apparent by considering a processes of crystallographic slip across $`K`$. To this end, imagine that a pair of representative atoms in $`K`$ undergo a relative sliding displacement of magnitude $`b`$ across a slip plane of the crystal. For this deformation one has

$$ϵ\left(K\right)=\frac{b}{h\left(K\right)}$$
(33)
It follows, therefore, that the adaption criterion (32) ensures that the mesh size $`h\left(K\right)<b`$ under the conditions just described. Thus, the adaption criterion is designed so that full atomistic resolution is attained when a simplex slips by a full Burgers vector. Interestingly, we have found that the dislocation patterns predicted in calculations are sensitive to the choice of tolerance, and generally $`TOL`$ must be chosen much smaller than $`1`$ in order not to inhibit dislocation nucleation.

## 3 Details of the computer implementation

In order for the quasicontinuum method to perform optimally, some implementational issues require careful attention. Several potential performance bottlenecks are discussed in this section.

### 3.1 Site-element mapping

The numerical solution of the quasicontinuum equilibrium equations (16), requires repeated computation of atomic coordinates $`𝒒_h\left(𝒍\right),𝒍`$ which depend, through the interpolation constraints (5), on the coordinates of representative atoms $`𝒒_h\left(𝒍_h\right),𝒍_h_h`$. It is therefore important to devise an efficient algorithm for locating the simplex $`K𝒯_h`$ which contains a given site $`𝒍`$ of the crystal.

The simplex $`K𝒯_h`$ containing site $`𝒍`$ may be found from elementary geometry and an exhaustive search over $`𝒯_h`$. Such method carries, however, substantial computational cost, due to a large number of required floating-point operations. Our approach is based on the fact that once the pair $`\{𝒍,K\}`$ has been established, it remains valid for the lifespan of $`𝒯_h`$. Thus, the following two stage strategy may be adopted:

1. A general search routine, which for a given lattice site $`𝒍`$ and triangulation $`𝒯_h`$ returns $`K𝒯_h`$ containing $`𝒍`$.
2. A look-up table (cache) that stores already associated pairs $`\{𝒍,K\}`$.

At first, the look-up table contains no $`\{𝒍,K\}`$ pairs and most of the inquiries results in the execution of the general search routine. As the force calculation proceeds, pairs are continuously inserted into the table, and fewer calls to the general routine are needed. At the end of the first force calculation, most of the pairs can be located in the cache, where they remain as long as $`𝒯_h`$ stays unchanged.

Our implementation of the look-up table is based on hashing, which is the method of referencing records in a table by doing arithmetic transformations of keys into table addresses. Any hashing algorithm requires two design decisions to be made:

1. A hash function $`h\left(k\right)`$ taking a key $`k`$ as its argument and returning an index into a table must be chosen.
2. One must establish a strategy for dealing with cases of two distinct keys $`k_1`$ and $`k_2`$ for which the resulting values of the hash function are the same, i. e. $`h\left(k_1\right)=h\left(k_2\right)`$ (collision resolution).

In his monograph, Knuth  provides an excellent description of hashing, including some possible choices of hash functions. For a simple three-dimensional Bravais lattice (1), the lattice coordinates of a site $`𝒍=\{l^1,l^2,l^3\}Z^3`$ may be considered as a key. The hash function acting on this key produces an index into a table of simplices where the $`K`$, containing $`𝒍`$, can be found. In cases in which a key contains more than one word (each of $`l^i`$ can be considered as a separate word) Knuth  suggests using the following hash function

$$h\left(𝒍\right)=\left[h_1\left(l^1\right)+h_2\left(l^2\right)+h_3\left(l^3\right)\right]modM.$$
(34)
where $`h_i\left(k\right)`$ denotes a hash function for $`l^i`$ and M is an integer parameter. The function $`amodp`$ computes the remainder of $`a/p`$. In general, the optimal choice of $`h_i`$ depends on $`𝒯_h`$ and cannot easily be established. We have found, however, that the particular hash function

$$h_i\left(l^i\right)=\text{lshift}(l^i,2^{i1}),i=1,2,3.$$
(35)
works well for a wide range of triangulations. The meaning of $`\text{lshift}(l,m)`$ is simply “shift all bits in $`l`$ left $`m`$ positions”. The operations involved in the calculation of $`h_i`$ are elementary, and the computational cost of each such calculation is small. The integer parameter $`M`$ in (34) may be chosen arbitrarily, but Knuth  points out that, when it is taken to be some power of $`2`$, i. e., $`M=2^p`$, $`amodM`$ is equivalent to masking the low $`p`$-bits from $`a`$ (on most of currently available computer hardware). The application of (34) to a lattice site $`𝒍`$ results in an index into the hash table from the interval $`[0,2^p)`$. Accordingly, the size of the table is limited to $`2^p`$, which renders collisions likely. The collisions are handled with the help of additional short tables which store $`\{𝒍,K\}`$ pairs for all $`𝒍`$ such that the $`h\left(𝒍\right)`$ produces the same value. Therefore, the search through the cache for a simplex $`K`$ containing site $`𝒍`$ becomes essentially a two stage process, in which the application of the hash function is followed by a linear search. Since, the number of lattice sites having the same index is small, and the additional cost due to the linear search phase is not significant.

### 3.2 Computation of cluster weights

As explained earlier, the cluster weights $`n_h\left(𝒍_h\right),𝒍_h`$ are obtained requiring that the shape functions are summed exactly by the cluster summation rule (17). The application of the summation rule to shape functions $`\phi _h\left(𝒍|𝒍_h\right)`$ leads to the system of $`N_h`$ linear algebraic equations

$$\underset{𝒍_h^{}_h}{}A\left(𝒍_h|𝒍_h^{}\right)n\left(𝒍_h^{}\right)=b\left(𝒍_h\right),𝒍_h_h,$$
(36)
to be solved for $`n_h\left(𝒍_h\right),𝒍_h`$. Here the matrix $`A\left(𝒍_h|𝒍_h^{}\right),𝒍_h^{},𝒍_h_h`$ is

$$A\left(𝒍_h|𝒍_h^{}\right)=\underset{𝒍^{}𝒞\left(𝒍_h^{}\right)}{}\phi _h\left(𝒍|𝒍_h\right).$$
(37)
and $`𝒞\left(𝒍_h^{}\right)`$ denotes the cluster of lattice sites centered at $`𝒍_h^{}`$. The array $`b\left(𝒍_h\right),𝒍_h`$ is simply

$$b\left(𝒍_h\right)=\underset{𝒍}{}\phi _h\left(𝒍|𝒍_h\right).$$
(38)
The calculation of the array $`b`$ requires the evaluation of the shape functions at all lattice sites within a crystal. While this may be regarded as acceptable for small crystals, it becomes prohibitively expensive when applied to large samples. However, as we have already shown, in the continuum limit summation may be replaced by integration, in which case simplex $`K`$ contributes $`N\left(K\right)/\left(d+1\right)`$ to each of its vertices. $`N\left(K\right)=\left(N/V\right)\left|K\right|`$ is the approximate number of sites within $`K`$ and $`N/V`$ is the atom density in the undeformed configuration of the crystal. Specifically, one may introduce a cutoff $`N_c`$ and restrict the direct calculation of the sum in (38) to simplices $`K`$ for which $`N\left(K\right)N_c`$.

It is apparent from (37) that $`A`$ has a sparse structure, which suggests the use of specialized solvers for sparse linear systems. However, a simpler alternative route is to resort to lumping in order to replace $`A`$ by a diagonal matrix (see, e. g., ). A lumping technique which is widely used in finite elements is the row-sum technique, which gives the diagonal entries of the lumped matrix as

$$A\left(𝒍_h|𝒍_h\right)=\underset{𝒍_h^{}_h}{}\underset{𝒍^{}𝒞\left(𝒍_h^{}\right)}{}\phi \left(𝒍|𝒍_h\right),𝒍_h_h,$$
(39)
with all other entries in $`A`$ set to zero. Once the matrix $`A`$ is lumped, the solution of (37) is trivial and gives

$$n\left(𝒍_h\right)=\frac{b\left(𝒍_h\right)}{A\left(𝒍_h|𝒍_h\right)},𝒍_h_h.$$
(40)
In the atomistic limit, the values of cluster weights computed from (40) become identical with obtained exactly.

## 4 Numerical tests

The accuracy of the quasicontinuum method is largely determined by three factors:

* The value of the cluster-size parameter $`r`$, which controls the accuracy of the cluster summation rules.
* The approximations introduced in the computation of the cluster weights, namely: the approximate calculation of vector $`b`$, controlled by the cuttoff parameter $`N_c`$; and the lumping procedure for constructing the diagonal matrix $`A`$.
* The value of parameter $`TOL`$ in (32), which controls the process of representative atom insertion.

In this section we present an analysis of the influence of each of these factors on the accuracy of the quasicontinuum method. We measure the quality of the solution in energy terms. Specifically, we identify the error in the approximate solution $`𝒒_h`$ with

$$e=\mathrm{\Phi }\left(𝒒_h\right)\mathrm{\Phi }\left(𝒒\right)$$
(41)
where $`𝒒`$ is the solution of the full system obtained, e. g., by relaxing $`𝒒_h`$. Since the potential energy decreases during this relaxation, the energy error is greater or equal to zero. For a well-supported harmonic crystal the energy error defined in (41) defines a proper norm of the error lattice function $`𝒒_h𝒒`$.

### 4.1 Test problem definition

We take *nanoindentation* as a convenient test problem for assessing the performance of the method. A salient feature of nanoindentation is the presence of a highly nonuniform state of deformation resulting in a sharp mesh gradation away from the indentor. This feature effectively tests the adaptivity of the method.

The test sample used in simulations is a cube of an fcc nearest-neighbor Lennard-Jones crystal containing $`64\times 64\times 64`$ fcc unit cells, or a total of $`\mathrm{1\hspace{0.17em}073\hspace{0.17em}345}`$ atoms. The size of the cube is limited by the need to compute the solution of the full problem by direct atomistic simulation in order to evaluate the error. The Lennard-Jones potential is

$$\varphi \left(r\right)=4\epsilon \left[\left(\frac{\sigma }{r}\right)^{12}2\left(\frac{\sigma }{r}\right)^6\right],$$
(42)
where $`\sigma `$ and $`\epsilon `$ are parameters which set the length and energy scales, respectively. The choice of the Lennard-Jones potential is motivated by the desire to eliminate surface effects which would otherwise compound the interpretation of results.

The surfaces of the sample are aligned with the cube directions, Fig. 6. Atoms in the bottom surface are constrained to remain at their initial positions throughout the test, whereas in all side surfaces the atoms are allowed to move in the $`z`$-direction only. No displacement constraints are introduced on the top surface of the cube.

In calculations we use a model of a spherical indentor proposed by Kelchner et al. . In this model, the indentor is regarded as an additional external potential $`\mathrm{\Phi }^{\mathrm{ext}}`$ interacting with atoms in the substrate. The potential has the particular functional form

$$\mathrm{\Phi }^{\mathrm{ext}}\left(r\right)=AH\left(Rr\right)\left(Rr\right)^3,$$
(43)
where $`R`$ is the radius of the indentor, $`r`$ denotes distance between a site and the center of the indentor, $`A`$ is a force constant and $`H\left(r\right)`$ is the step function. In calculations we adopt the following values of the parameters: $`R=100\left[\sigma \right]`$ and $`A=2000\left[\epsilon /\sigma ^3\right]`$.

The initial triangulation of the cube is specifically tailored to the nanoindentation geometry, Fig. 7. Thus, in the region of the crystal located directly underneath the indentor, full atomistic resolution is introduced from the outset. Away from this region, the triangulation becomes gradually coarser. The resulting number of representative atoms in the initial mesh is $`888`$, a significant reduction from the total number of atoms ($`\mathrm{1\hspace{0.17em}073\hspace{0.17em}345}`$) in the sample. A cross section of the initial mesh through the center of the cube with a plane $`y=const`$ is show in Fig. 7. All solutions are computed using a nonlinear version of the Conjugate Gradient method.

### 4.2 Effect of the cluster size

The cluster radius $`r`$ controls the accuracy of the cluster summation rules. In order to isolate the effect of the cluster size from other factors, in the present tests the summation weights are computed exactly and the mesh is kept fixed. The indentor is pushed through a distance $`0.1\left[\sigma \right]`$ into the crystal along the z-axis. This value ensures that the material immediately under the indentor becomes highly deformed and is pushed well into the nonlinear regime. However, the induced deformation is not sufficient for any defects to arise. This greatly simplifies the interpretation of results and the error analysis, since under the stated conditions the full atomistic solution $`𝒒`$ to compare with is well-defined.

The computed dependence of the error $`e`$ on the cluster size $`r`$ is shown in Fig. 8. It is evident from this plot that the inclusion of a single shell of neighbors in the clusters suffices to eliminate the rank-deficiency of the node-based summation rule. As suggested by our previous analysis, the error is greatest for the smallest cluster size and decreases rapidly as the cluster size increases. For $`r`$ above $`\sqrt{8}\left[\sigma \right]`$, corresponding to the radius of the eighth shell of neighbors, the solution becomes insensitive to the cluster size. Thus, the use of relative small clusters results in high accuracy while preserving the computational efficiency of the method.

### 4.3 Effect of the lumping procedure

The error in the calculation of the cluster weights originates from two sources: the approximate computation of the array $`b`$ in (38), controlled by the value of parameter $`N_c`$; and the lumping of matrix $`A`$ in (37). In order to appraise this error, we repeat the calculations described in the preceding section with the summation weights computed approximately. The cutoff $`N_c`$ is set to $`2000`$, i. e. the contributions to the array $`b`$ from simplices which contain fewer than $`2000`$ atoms are computed explicitly. The dependence of the error $`e`$ on the cluster size parameter $`r`$ for both the approximate (squares) and exact (circles) weights is shown in Fig. 9. Small differences do arise for very small clusters, and virtually disappear for large clusters. In view of these results we may conclude that the use of approximate summation weights does not result in a significant loss of accuracy.

### 4.4 Convergence properties of the quasicontinuum method

We proceed to assess numerically the rate at which the quasicontinuum solution converges towards the fully atomistic solution. The test case is as in the preceding sections. In all cases, the cluster size $`r`$ is set to $`2\sqrt{2}\left[\sigma \right]`$ and the cluster summation weights are computed approximately with a cutoff parameter $`N_c=2000`$. The coarsest mesh is shown in Fig. 10a. Subsequent finer meshes are constructed by regular refinement, Fig. 10b, resulting in increasing numbers of representative atoms $`N_h`$, with the exception of the fully resolved case, $`h=\sqrt{2}\left[\sigma \right]`$, in which the representative atoms are placed on all the crystal lattice sites. In these calculations, the crystal is loaded simply by imparting a downward displacement $`0.1\left[\sigma \right]`$ to the central atom on the surface labeled $`A`$ in Fig. 6.

Fig. 11 shows the variation of the energy error with the number $`N_h`$ of representative atoms. Three distinct regimes are evident in this figure. For small values of $`N_h`$, the error behaves as $`eO\left(N_{h}^{}{}_{}{}^{\alpha }\right)`$, with a convergence rate $`\alpha 0.39`$. At $`N_h=4913`$, corresponding to $`h=4\sqrt{2}\left[\sigma \right]`$, the summation clusters begin to overlap and the rate of convergence increases to $`\alpha 0.55`$. The final drop of the error occurs when full atomistic resolution is attained.

These results demonstrate the convergence characteristics of the quasicontinuum method in the extreme case of an applied point load. When full lattice sums are performed the calculated convergence rate of $`0.55`$ approaches the optimal rate of convergence of linear finite elements applied to linear elasticity, which is $`2/3`$. The introduction of cluster summation rules results in a degradation of the convergence rate. It should be carefully noted, however, that some of these conclusions may well depend on the loading geometry. As is well known, the linear elasticity solution corresponding to a point load diverges under the point of application of the load and does not possess finite energy. The finiteness of the energy of atomistic solution owes entirely to the discreteness of the atomic lattice. In this sense, the test case considered here is, therefore, a worse case, and it is possible that the converge of the quasicontinuum method is more robust for smooth loading.

### 4.5 Effect of the remeshing-indicator tolerance

The lack of convexity of the energy functional for crystalline materials allows for lattice defects to develop. Among these defects, dislocations play a crucial role, as the carriers of the plastic deformation. One of the appealing features of the quasicontinuum method is its ability to follow the nucleation and motion of dislocations inside the crystal. This is achieved by adaptively providing full atomistic resolution in the highly deformed regions of the crystal. In the nanoindentation test, dislocations nucleate under the indentor upon the attainment of a critical load and subsequently propagate into the crystal.

The insertion of new representative atoms, leading to mesh refinement, is controlled by criterion (32) and, specifically, by the tolerance $`TOL`$. A small value of $`TOL`$ promotes refinement, whereas a large value of $`TOL`$ inhibits refinement. Most importantly, an excessively large value of $`TOL`$ may have the undesirable effect of inhibiting the nucleation of dislocations under the indentor altogether. On the other hand, an unduly small value of $`TOL`$ results in excessive refinement and a prohibitive computational burden.

In this section we explore this trade-off by way of numerical testing, with a view to bracketing the optimal value of the tolerance parameter $`TOL`$. To this end, we fix the cluster size to $`r=\sqrt{2}\left[\sigma \right]`$, and we use the approximate cluster summation weights with a cutoff $`N_c=2000`$. The indentor is driven into the crystal at increments of $`0.1\left[\sigma \right]`$ up to a maximum indentation depth of $`2.0\left[\sigma \right]`$. The first loading step does not result in the nucleation of dislocations, even in the fully atomistic model, and is identical to the calculations described in the preceding sections. By contrast, at the indentation depth of $`2.0\left[\sigma \right]`$ the presence of dislocation structures under the indentor is clearly evident in the fully atomistic model, Fig. 12.

The energy error $`e`$ is plotted in Fig. 13 as a function of indentation depth $`\delta `$ for different values of the tolerance parameter $`TOL`$. Initially, the error is insensitive to $`TOL`$, and all curves ostensibly follow the same linear growth pattern. However, the error curves diverge after a penetration depth $`\delta =0.5\left[\sigma \right]`$. For large values of $`TOL`$, the error continues to grow owing to the inability of the model to resolve the emerging dislocation structures. Indeed, the nucleation of dislocations during the test is entirely inhibited by all but the smallest values of the tolerance, $`TOL=0.004`$. For small values of the tolerance, $`TOL=0.005`$ and $`0.004`$, the energy error reaches a plateau at roughly $`1\%`$ of the total energy, and remains constant or decreases thereafter.

The final mesh at $`\delta =2.0\left[\sigma \right]`$ corresponding to a tolerance of $`TOL=0.004`$ is shown in Fig. 14. The total number of representative atoms in this configuration is $`N_h=141,469`$, which corresponds approximately to $`1/8`$ of the total number $`N`$ of atoms in the sample. It should be stressed, however, that since all of new representative atoms are inserted in the vicinity of the indentor, the ratio $`N_h/N`$ can be made arbitrarily small by considering increasingly large samples.

The dislocation structures predicted by the quasicontinuum method for a tolerance $`TOL=0.004`$, Fig. 15, should be contrasted with those predicted by the fully atomistic model, Fig. 12. The figures display the energetic atoms (red) underneath the surface of the crystal (blue). Evidently, these structures differ in detail. This is expected, owing to the massive lack of uniqueness which characterizes the problem. Thus, as dislocation are nucleated and the lack of convexity of the energy comes into play, numerous deformation paths become available to the crystal. Many of these deformation paths differ only slightly in their energy content and are, therefore, indistinguishable at the macroscale, e. g., in terms of the corresponding force-depth of penetration curves. This is specially so when the crystal is loaded along a direction of high symmetry, which results in the possible activation of a large number of competing slip systems. Under these conditions, a pointwise comparison of solutions becomes mathematically meaningless, and the sole meaningful criterion to measure the quality of a solution is its energy contents. By this criterion, the quasicontinuum and fully atomistic solutions are of comparable quality, as their energies differ by less than $`1\%`$.

This contention is confirmed by Fig. 16, which compares the force on the indentor $`f`$ as a function of penetration depth $`\delta `$ predicted by the fully atomistic and quasicontinuum ($`TOL=0.004`$) models. Indeed, the two curves are ostensibly indistinguishable even beyond the onset of dislocation nucleation. On this basis, in practice choices of the tolerance $`TOL`$ in the range of $`0.001`$$`0.004`$ would appear to strike an adequate balance between accuracy and performance demands.

## 5 Summary and Discussion

We have developed a streamlined and fully three-dimensional version of the quasicontinuum method of Tadmor et al.  and we have presented a numerical analysis of its accuracy and convergence characteristics. As a new addition to the theory, we have formulated a new class of summation rules in which the lattice function being summed is sampled over clusters of atoms. The size of these clusters is an adjustable parameter and controls the accuracy of the summation rule. We have also presented an efficient method for computing the requisite summation weights in linear time. Beyond its usefulness as a numerical scheme for approximating lattice sums, the cluster approach provides a systematic means of sampling the behavior of small representative crystallites, and thus opens a possible avenue for incorporating additional physics to quasicontinuum models such as diffusion of solute atoms and thermal lattice vibrations.

We have presented a suite of numerical tests which demonstrate the accuracy and performance of the method. As expected, the accuracy of cluster summation rules increases with cluster size. Furthermore, our numerical tests suggest that the addition of a single shell of neighbors suffices to stabilize the rank-deficiency which afflicts node-based summation rules. The computed convergence rate of the method in problems involving the application of point loads to crystals is close to linear when the lattice sums are performed exactly, and decreases somewhat when the sums are approximated using a cluster summation rule. It is worth noting that, contrary to continuum problems for which a well-developed approximation theory exists, a similar approximation theory for lattice problems appears to be missing at present, even for linear problems. The development of rigorous error bounds for finite-element approximations to lattice problems is a clear worthwhile direction of future research.

Perhaps the most interesting issue among those addressed in this paper concerns the ability of the quasicontinuum method to simulate microstructural evolution, which owes largely to mesh adaptivity. The same lack of convexity which allows for defects and microstructures to arise in the first place renders solutions massively nonunique. For a given loading or prescribed deformation, it is in general possible to find multiple equilibrium states of a crystal with defects possessing equal or nearly equal energies. These states are indistinguishable from a macroscopic point of view, e. g., in the sense of yielding identical force-displacement curves. Whether the solution follows one deformation path or another depends sensitively on small perturbations of the system, including details of the mesh design. Under these conditions, the pointwise comparison of solutions is not particularly meaningful. Our numerical tests suggest that, with sufficient mesh adaptivity, the quasicontinuum method is capable of simulating evolving microstructures comparable, in energetic terms, to those obtained from a full atomistic calculation.

## Acknowledgments

The support of the Department of Energy through Caltech’s ASCI/ASAP Center for the Simulation of the Dynamic Behavior of Solids is gratefully acknowledged.