text:
# A Bernstein-Bézier Sufficient Condition for Invertibility of Polynomial Mapping Functions (Draft)

## 1 Invertibility of polynomial mapping functions

In finite element analysis, it is common to subdivide the domain into elements that are images of a reference domain under polynomial functions. This approach gives rise to the popular isoparametric elements . The reference domain is the unit square $`I^2=\{(\xi ,\eta ):0\xi ,\eta 1\}`$ or triangle $`\mathrm{\Delta }^2=\{(\xi ,\eta ):0\xi ,\eta ;\xi +\eta 1\}`$ in $`𝐑^2`$ or the unit cube $`I^3=\{(\xi ,\eta ,\zeta ):0\xi ,\eta ,\zeta 1\}`$ or tetrahedron $`\mathrm{\Delta }^3=\{(\xi ,\eta ,\zeta ):0\xi ,\eta ,\zeta ;\xi +\eta +\zeta 1\}`$ in $`𝐑^3`$.

Polynomials defined on $`\mathrm{\Delta }^d`$, $`d=2,3`$, generally include monomial terms up to degree $`p`$ for some $`p>0`$. On the other hand, for the unit cube $`I^d`$, the monomial terms are generally up to degree $`p`$ individually in each coordinate. Therefore, for the rest of the paper we use $`p`$ to denote the maximum total degree of polynomials on $`\mathrm{\Delta }^2`$ or $`\mathrm{\Delta }^3`$, and the maximum degree in individual coordinates of polynomials defined on $`I^2`$ or $`I^3`$.

Let $`F:U𝐑^d`$ be a polynomial function, where $`U`$ is either $`I^d`$ or $`\mathrm{\Delta }^d`$ and $`d`$ is 2 or 3. Suppose $`F`$ is injective on $`U`$, implying that it has an inverse $`G`$ defined on $`F(U)`$ such that $`G(F(𝐮))=𝐮`$. If such a $`G`$ exists and is smooth, then we say that $`F`$ is invertible. Invertibility is important for correctness of a finite element mesh. This definition of invertibility is sometimes called “global invertibility” in order to distinguish it from “local invertibility.” The function $`F`$ is locally invertible if its derivative, denoted as $`J`$, is nonsingular on the entire reference element $`U`$. Local invertibility is necessary for global invertibility but is not sufficient. A function can be locally but not globally invertible if $`F(U)`$ “wraps around” (e.g., consider $`F(r,\theta )=r(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$ for $`(r,\theta )[1,2]\times [0,4\pi ]`$). Additional sufficient conditions for global invertibility in the general setting are considered by Ivanenko . Beyond local invertibility, one would also like an upper bound on the “condition number”

$$cond(F)=\underset{𝐮U}{\mathrm{max}}J(𝐮)\underset{𝐮U}{\mathrm{max}}J(𝐮)^1$$
(1)
of $`F`$. Even more strongly, one would like bounds on the higher derivatives of $`F`$ in terms of $`cond(F)`$ to satisfy the Ciarlet-Raviart  conditions on convergence of order-$`p`$ finite element approximation.

Unfortunately, even the first problem, namely, determining nonsingularity of the Jacobian, is a difficult problem. We know of simple necessary and sufficient conditions only for the following special cases:

1. linear polynomials, that is, $`p=1`$, defined on $`\mathrm{\Delta }^d`$ of all dimensions,
2. quadratic polynomials on $`\mathrm{\Delta }^2`$ in the case $`d=2`$, in the special case that $`F`$ is linear on two of the three edges of the reference triangle , and
3. bilinear elements (i.e., $`p=1`$) defined on $`I^2`$ .

For all other cases, we know of only separate necessary and sufficient conditions, and this paper also proposes only a sufficient condition. To our knowledge, the only necessary condition for the general case proposed in the literature is that $`det(J)`$ have the same sign (strictly positive or strictly negative) at some finite list of test points. In engineering applications, it is common to test invertibility at the Gauss points used for quadrature on the element . This test is known not to be sufficient; see further remarks on this point below.

Sufficient conditions have been proposed for a few settings. For example Ushakova  presents a sufficient condition for the case $`p=1`$ on domain $`I^3`$. Her sufficient condition is interpreted as requiring that certain tetrahedra formed by choosing subsets of four of the eight vertices of $`F(I^3)`$ have positive volume. Sufficient conditions for quadratic triangles and tetrahedra have been proposed by Salem, Canann and Saigal .

In this work, we present a sufficient condition to ensure both local and global invertibility of $`F`$ for any of the four domains $`I^d`$, $`\mathrm{\Delta }^d`$, $`d=2,3`$ and for any polynomial degree $`p`$. Our condition appears to be weaker (i.e., not able to certify invertibility in more cases) than existing sufficient tests for specific reference domains and values of $`p`$ but is considerably more general.

Our condition is based on writing $`J`$ in Bernstein-Bézier form and then considering the convex hull of the derivatives at the control points. In the next section we describe the new sufficient condition and establish that it is sufficient for local invertibility. In Section 3 we provide an equivalent characterization of the sufficient condition that is amenable to efficient testing. In Section 4 we use the equivalent characterization to prove that the condition also implies global invertibility. Our condition has a desirable property that we term “affinity,” which we define in the last section.

## 2 Bernstein-Bézier form

Bernstein-Bézier (BB) form is a popular way to write polynomials in computer-aided geometric design . A univariate polynomial of degree $`p`$ in BB form would be written:

$$F(\xi )=\underset{i=0}{\overset{p}{}}f_i(1\xi )^{pi}\xi ^i\frac{p!}{i!(pi)!}$$
(2)
and has natural parametric domain $`\xi [0,1]`$. A bivariate polynomial with maximum degree $`p`$ individually in $`\xi ,\eta `$ is written in the form:

$$F(\xi ,\eta )=\underset{i=0}{\overset{p}{}}\underset{j=0}{\overset{p}{}}𝐟_{i,j}(1\xi )^{pi}\xi ^i(1\eta )^{pj}\eta ^j\frac{p!p!}{i!(pi)!j!(pj)!}$$
(3)
and has as its natural parametric domain $`(\xi ,\eta )I^2`$. A trivariate polynomial with maximum degree $`p`$ individually is written

$$F(\xi ,\eta ,\zeta )=\underset{i=0}{\overset{p}{}}\underset{j=0}{\overset{p}{}}\underset{k=0}{\overset{p}{}}𝐟_{i,j,k}(1\xi )^{pi}\xi ^i(1\eta )^{pj}\eta ^j(1\zeta )^{pk}\zeta ^k\frac{p!p!p!}{i!(pi)!j!(pj)!k!(pk)!},$$
(4)
and has as its natural parametric domain $`(\xi ,\eta ,\zeta )I^3`$.

A bivariate polynomial with total degree at most $`p`$ is written in the BB form

$$F(\xi ,\eta )=\underset{i=0}{\overset{p}{}}\underset{j=0}{\overset{pi}{}}𝐟_{i,j}\xi ^i\eta ^j(1\xi \eta )^{pij}\frac{p!}{i!j!(pij)!}$$
(5)
and has natural parametric domain $`(\xi ,\eta )\mathrm{\Delta }^2`$.

Finally, a trivariate polynomial with total degree at most $`p`$ is written in the BB form

$$F(\xi ,\eta ,\zeta )=\underset{i=0}{\overset{p}{}}\underset{j=0}{\overset{pi}{}}\underset{k=0}{\overset{pij}{}}𝐟_{i,j,k}\xi ^i\eta ^j\zeta ^k(1\xi \eta \zeta )^{pijk}\frac{p!}{i!j!k!(pijk)!}$$
(6)
and has natural parametric domain $`(\xi ,\eta ,\zeta )\mathrm{\Delta }^3`$.

In all five cases, the vectors $`f_i`$ (in 1D), $`𝐟_{i,j}𝐑^2`$ (in 2D) or $`𝐟_{i,j,k}𝐑^3`$ (in 3D) are called the control points. A fundamental theorem about BB form is:

###### Theorem 1

Let $`U`$ be the natural parametric domain of BB form for the five cases listed above. Then $`F(U)`$ is contained in the convex hull of the control points.

Farin  proves this as a consequence of the deCasteljau algorithm for evaluating $`F`$, but here is a sketch of a more direct proof. One observes that on the natural parametric domain, all the factors in the summations like $`\xi `$, $`(1\xi \eta )`$, etc., are nonnegative. Furthermore, one observes that the value of all the sums, if the control points are excluded, is exactly 1. For example, consider $`(\text{2})`$ without control points:

$$\underset{i=0}{\overset{p}{}}(1\xi )^{pi}\xi ^i\frac{p!}{i!(pi)!}$$

This sum is identically 1, as seen by considering a binomial expansion of $`(\xi +(1\xi ))^p`$. Similar argument apply to $`(\text{3})`$$`(\text{6})`$. Thus, in the natural parametric domains, the above summations may be regarded as weighted averages of the control points, proving the theorem.

The next feature of BB form is that once a function is in BB form, the derivatives can be easily put into BB form. In all cases, the derivative is a BB expansion of one degree lower, multiplied by $`p`$, and with control points that are finite differences of the control points for $`F`$ in the direction of the variable being differentiated. Thus, for example, in the case of $`(\text{5})`$, we have

$$\frac{F}{\xi }=\underset{i=0}{\overset{p1}{}}\underset{j=0}{\overset{pi1}{}}p(𝐟_{i+1,j}𝐟_{i,j})\xi ^i\eta ^j(1\xi \eta )^{pij1}\frac{(p1)!}{i!j!(pij1)!}$$

and

$$\frac{F}{\eta }=\underset{i=0}{\overset{p1}{}}\underset{j=0}{\overset{pi1}{}}p(𝐟_{i,j+1}𝐟_{i,j})\xi ^i\eta ^j(1\xi \eta )^{pij1}\frac{(p1)!}{i!j!(pij1)!}.$$

In this case, the list of $`p(p+1)/2`$ vectors of the form $`p(𝐟_{i+1,j}𝐟_{i,j})`$ are control points of $`F/\xi `$, and analogously $`p(𝐟_{i,j+1}𝐟_{i,j})`$ are control points for $`F/\eta `$. Let us denote these two lists of control points $`G_\xi `$ and $`G_\eta `$ respectively. Similar expressions hold for the other four BB forms described above. In 3D there is a third list $`G_\zeta `$. We now state our condition:

###### Condition 1

In the case $`d=2`$, the matrix $`[𝐮,𝐯]`$ is invertible for every $`𝐮hull(G_\xi )`$ and every $`𝐯hull(G_\eta )`$. In the case $`d=3`$, the matrix $`[𝐮,𝐯,𝐰]`$ is invertible for every $`𝐮hull(G_\xi )`$, $`𝐯hull(G_\eta )`$, $`𝐰hull(G_\xi )`$.

In this condition, “hull” denotes the convex hull. The main result of this section is as follows.

###### Theorem 2

If Condition 1 holds, then the matrix $`J`$ is invertible on the entire reference element.

The proof of this theorem follows from the arguments in the previous paragraphs: Condition 1 is sufficient since the actual Jacobians that occur on the domain have their first columns chosen from $`hull(G_\xi )`$, etc.

## 3 A computational characterization of the sufficient condition

Condition 1 specifies our sufficient condition in somewhat nonconstructive terms. In this section we provide a computational means to verify the sufficient condition.

###### Theorem 3

In the case of $`𝐑^2`$, Condition 1 is equivalent to the following condition:

* there exists a vector $`𝐡𝐑^2`$ such that for all $`𝐟G_\xi G_\eta `$, $`𝐡^T𝐟>0`$, and
* there exists a vector $`\overline{𝐡}𝐑^2`$ such that for all $`𝐟G_\xi `$, $`\overline{𝐡}^T𝐟>0`$ and for all $`𝐟G_\eta `$, $`\overline{𝐡}^T𝐟<0`$.

* Proof. Let $`𝐮,𝐯`$ be arbitrary in $`hull(G_\xi )`$ and $`hull(G_\eta )`$ respectively. This is equivalent to saying there exist nonnegative $`\alpha _1,\mathrm{},\alpha _n`$ summing to 1 such that $`𝐮=\alpha _1𝐟_1+\mathrm{}+\alpha _n𝐟_n`$, where $`𝐟_1,\mathrm{},𝐟_n`$ is an enumeration of $`G_\xi `$. Similarly, there exist nonnegative $`\beta _1,\mathrm{},\beta _m`$ summing to 1 such that $`𝐯=\beta _1𝐠_1+\mathrm{}+\beta _m𝐠_m`$, where $`𝐠_1,\mathrm{},𝐠_m`$ is an enumeration of $`G_\eta `$. The condition that $`[𝐮,𝐯]`$ is invertible is equivalent to saying that $`𝐮,𝐯`$ are independent, i.e., that there do not exist $`\delta ,\gamma `$, at least one nonzero such that $`\delta 𝐮+\gamma 𝐯=\mathrm{𝟎}`$. Suppose they are dependent, and let $`\delta ,\gamma `$ be the two coefficients of dependence. Without loss of generality, $`\delta 0`$. There are now two cases: either $`\gamma 0`$ or $`\gamma <0`$. If $`\gamma 0`$, define $`\overline{\alpha }_1,\mathrm{},\overline{\alpha }_n`$ to be $`\delta \alpha _1,\mathrm{},\delta \alpha _n`$ and define $`\overline{\beta }_1,\mathrm{},\overline{\beta }_m`$ to be $`\gamma \beta _1,\mathrm{},\gamma \beta _n`$. Then the condition that the $`\alpha `$’s are nonnegative and the assumption that $`\delta 0`$ is equivalent to the hypothesis that $`\overline{\alpha }_1,\mathrm{}\overline{\alpha }_n`$ are nonnegative (with no restriction on their sum). Similarly, the condition on the $`\overline{\beta }_i`$’s is that they are nonnegative.

Thus, in the case that $`\gamma 0`$, a dependence is equivalent to the existence of nonnegative $`\overline{\alpha }_i`$’s and $`\overline{\beta }_i`$’s, not all zeros, such that

$$\underset{i=1}{\overset{n}{}}\overline{\alpha }_i𝐟_i+\underset{i=1}{\overset{m}{}}\overline{\beta }_i𝐠_i=\mathrm{𝟎}.$$

This problem is solved by linear programming. By Farkas’ lemma , there is a nonnegative solution to this problem, with not all the coefficients zero, if and only if there does not exist a vector $`𝐡`$ such that $`𝐡^T𝐟_i>0`$ and $`𝐡^T𝐠_i>0`$ for all $`𝐟_i`$’s and $`𝐠_i`$’s.

Now we turn to the case of dependence when $`\gamma <0`$. In this case, using the analogous argument, there is no dependence provided that there is a vector $`\overline{𝐡}`$ such that $`\overline{𝐡}^T𝐟_i>0`$ for all $`i`$ and $`\overline{𝐡}^T(𝐠_i)>0`$ for all $`i`$.

This shows that in the case of $`𝐑^2`$, Condition 1 can be tested by solving two linear programming problems over $`𝐑^2`$ to find $`𝐡`$ and $`\overline{𝐡}`$. In fact, there is a much simpler algorithm, namely, compute the args of all the points in $`G_\xi G_\eta `$. To determine whether $`𝐡`$ exists, one checks whether the max arg differs by less than $`\pi `$ from the min arg. (The branch cut for defining the arg function must be chosen outside the min-max range). A similar test can determine whether $`\overline{𝐡}`$ exists.

In the three-dimensional setting, the algorithm is not so simple but it is still linear time. Following the same approach as in the previous proof, we see that in 3D, Condition 1 is equivalent to the following conditions:

* there exists a vector $`𝐡_1𝐑^3`$ such that $`𝐡_1^T𝐟>0`$ for any $`𝐟G_\xi G_\eta G_\zeta `$,
* there exists a vector $`𝐡_2𝐑^3`$ such that $`𝐡_2^T𝐟>0`$ for any $`𝐟G_\xi G_\eta (G_\zeta )`$, where $`G_\zeta `$ means $`\{𝐟:𝐟G_\zeta \}`$,
* there exists a vector $`𝐡_3𝐑^3`$ such that $`𝐡_3^T𝐟>0`$ for any $`𝐟G_\xi (G_\eta )G_\zeta `$, and
* there exists a vector $`𝐡_4𝐑^3`$ such that $`𝐡_4^T𝐟>0`$ for any $`𝐟G_\xi (G_\eta )(G_\zeta )`$.

Each of these can be verified with linear programming, which is linear time in three dimensions . Alternatively, they can be verified in $`O(n\mathrm{log}n)`$ time using a convex hull algorithm .

## 4 Global invertibility

We are now in a position to prove that Condition 1 in fact implies global invertibility of the mapping function. To demonstrate global invertibility, we will show that $`F`$ is injective, i.e., for all $`𝐮,𝐯`$ in the reference domain, $`𝐮𝐯F(𝐮)F(𝐯)`$. Injectivity is sufficient for global invertibility: if an injective smooth function has a nonsingular derivative at all points of its domain, then the global inverse is also smooth by the inverse mapping theorem. The nonsingularity of the derivative is already established.

###### Theorem 4

If Condition 1 holds, then $`F`$ is globally invertible on the reference element.

* Proof. Let us start with the proof of injectivity in two dimensions. Choose $`𝐮,𝐯`$ in the reference element such that $`𝐮𝐯`$. Let $`𝐰=𝐯𝐮`$. Then by the definition of directional derivative,

$$F(𝐯)F(𝐮)=_0^1J(𝐮+t𝐰)𝐰𝑑t$$
(7)
Assume without loss of generality that $`w_10`$, where $`(w_1,w_2)`$ denotes the entries of $`𝐰`$. (The case $`w_1<0`$ is handled by exchanging the roles of $`𝐮,𝐯`$.) Now there are two cases: either $`w_20`$ or $`w_2<0`$. Suppose first that $`w_20`$. Because we assume Condition 1 and because we have established Theorem 3, we conclude that there exists a vector $`𝐡𝐑^2`$ such that $`𝐡^T𝐟>0`$ for all $`𝐟G_\xi G_\eta `$. By convexity, this implies $`𝐡^T𝐟>0`$ for all $`𝐟hull(G_\xi G_\eta )`$. This implies that $`𝐡^TJ(\xi ,\eta )`$ is a vector both of whose entries are positive for any $`\xi ,\eta `$ in the reference domain since both columns of $`J`$ are taken from $`hull(G_\xi G_\eta )`$. Since $`𝐰`$ has both entries nonnegative, and at least one entry of $`𝐰`$ is positive (because $`𝐮𝐯`$), this means $`𝐡^TJ(\xi ,\eta )𝐰>0`$ for all $`\xi ,\eta `$ in the reference domain. Thus, taking the inner product of both sides of $`(\text{7})`$ with $`𝐡`$ shows that $`𝐡^T(F(𝐯)F(𝐮))>0`$, and in particular, $`F(𝐯)F(𝐮)`$.

The second case is $`w_2<0`$. The argument is analogous, except we use $`\overline{𝐡}`$ instead of $`𝐡`$.

Finally, the extension to 3D follows the same lines. We assume $`w_10`$. Depending on the signs of entries $`w_2,w_3`$ of $`𝐰`$, we select one of $`𝐡_1,𝐡_2,𝐡_3,𝐡_4`$.

## 5 Discussion

We say that a sufficient condition for invertibility has the “affinity” property provided that it is always satisfied if the mapping function is a nondegenerate affine linear mapping, or a sufficiently small perturbation of a nondegenerate affine linear mapping.

Affinity is a desirable property since the affine linear mapping is obviously the easiest case for invertibility. Furthermore, if a sufficient condition for invertibility has the affinity property, then it can be turned into a necessary and sufficient condition if applied to each subcell of a sufficiently fine subdivision of the original domain. This is because an invertible mapping on a fine subdivision will behave like an affine mapping plus a small perturbation on each subcell, and the small perturbation tends to zero as the subdivision gets finer.

It is fairly easy to see that our test has the affinity property. Recall that the BB control points for a constant function are all identical. Thus, if $`F`$ is affine linear, then $`G_\xi `$, $`G_\eta `$ and $`G_\zeta `$ are all singleton sets. Assuming that $`F`$ is affine and invertible means that the condition must be satisfied. Furthermore, a sufficiently small perturbation of $`F`$ will give $`G_\xi `$, $`G_\eta `$, and $`G_\zeta `$ a positive radius, but the condition of separability by planes will still hold.

Ushakova’s condition also has the affinity property. We suggest that any proposed sufficient condition should possess the affinity property in order to be considered useful.

Ushakova’s sufficient condition for the $`p=1`$ and $`U=I^3`$ case is less restrictive (i.e., is able to validate more elements) than the condition proposed here.

We conclude with a few open questions raised by this work.

1. In the case $`p=1`$ and $`U=I^3`$, our condition is more restrictive than Ushakova’s. On the other hand, it is not obvious how to generalize her condition to higher degrees or to tetrahedral elements. It would be interesting to combine her techniques with the ones in this paper to come up with less restrictive conditions for higher degrees. Ushakova (private communication) found that experiments with randomly generated elements indicate that for $`p=1`$, $`U=I^3`$, the conditions stated here overly restrictive in the sense that most of the valid (invertible) elements among a randomly generated test set would not satisfy the sufficient condition proposed here.
2. In practical application of isoparametric mesh generation, one defines the polynomial function only on one face of the boundary of the element (namely, the face adjacent to a curved exterior boundary). The remaining coefficients defining $`F`$ can be determined by a formula such as the formula proposed by Lenoir , which has certain theoretical guarantees. How does the sufficient condition here specialize if we assume that interior degrees of freedom of $`F`$ are determined by Lenoir’s formula? Our very preliminary computational tests seem to indicate that using Lenoir’s formula seems to make the element more amenable to our sufficient condition (i.e., it seems to minimize the diameters of the sets $`G_\xi ,G_\eta ,G_\zeta `$ compared to other choices for interior degrees of freedom).
3. The test proposed here could be strengthened to obtain an upper bound on $`cond(F)`$ defined by $`(\text{1})`$ in terms of the numerical values of $`𝐡^T𝐟_i`$ as $`𝐟_i`$ ranges over the sets $`G_\xi ,G_\eta ,G_\zeta `$ and $`𝐡`$ is the vector defining one of the halfspaces in Theorem 3. This raises the possibility of a sufficient condition not only for invertibility but also for confirming the Ciarlet-Raviart conditions. Recall that the Ciarlet-Raviart conditions require that higher derivatives of $`F`$ are bounded in terms of $`cond(F)`$.

Note that other condition numbers besides $`cond(F)`$ may be useful in practice. For example, Branets and Garanzha  have developed a distortion measure related to structured grid generation. It would be interesting to get bounds on all of these condition numbers in terms of the the numerical values of $`𝐡^T𝐟_i`$.

## 6 Acknowledgments

The author is grateful to O. Ushakova, S. Ivanenko and V. Garanzha for helpful comments on this manuscript.


format: text

text:
# Sampling the diffusion paths of a neutral vacancy in Silicon with SIEST-A-RT

## I Introduction

Self-diffusion in homogeneous solids is a fundamental process of mass transport, in addition, it is responsible of the annealing of implantation damage, crystal to amorphous transition and the nucleation of extended defects. The native defects, such as vacancies and self-interstitial, have been identified as the prime entities controlling the self-diffusion, as well as the diffusion of impurities in semiconductor lattices. For example, in silicon, antimony impurities are vacancy diffusers and phosphorus are interstitial diffusers Ural et al. (1999). Therefore, self-diffusion in Si has been extensively investigated by experiments Lagna and Coffa (1999); Watkins (1992); Watkins and Corbett (1964); Dannefaer et al. (1986); Ural et al. (1999) and simulations Ghaisas (1991); Öğüt and Chelikowsky (1999); Jääskeläinen. et al. (2001); Dobson et al. (1989); Seong and Lewis (1996); Puska et al. (1998); Antonelli et al. (1998); Bar-Yam and Joannnopoulos (1984); Kumeda et al. (2001); Virkkunen et al. (1993); Goedecker et al. (2002); Marques et al. (2001); Kim et al. (2000); Wang et al. (1991); Tang et al. (1997); Leung et al. (1999); Estreicher (2000); Henkelman et al. (2002); Munro and Wales (1999); Puska (2000); Mercer et al. (1998); Probert and Payne (2003); Henkelman et al. (2000); Roberson and Estreicher (1994); Öğüt et al. (1997); Anderson et al. (1996); Clarck and Ackland (1997); Lento and Nieminen (2003).

In covalently bonded silicon, the neutral vacancy is the most stable charge state near room temperature Watkins (1992); Agullo-Lopez et al. (1988). While the diffusion mechanism and the barrier for silicon have been experimentally determined decades ago Dannefaer et al. (1986); Watkins and Corbett (1964), we still lack a detailed knowledge of the various activated mechanisms associated with the neutral vacancy. With standard ab initio molecular dynamics methods, the number of events generated is always insufficient to give reliable data for migration energies. Accordingly, most of the focus has been on formation energies, distortion and reorientation geometries Anderson et al. (1996). Few papers reported a calculated activation energy and a migration path that are in agreement with the experimental data Roberson and Estreicher (1994); Maroudas and Brown (1993). Until now, the only studies available for exploring the energy landscape were that of Munro and Wales Munro and Wales (1999) who examined the dynamics of a vacancy, among others, using a tight-binding approach, and the work of Kumeda et al. Kumeda et al. (2001) using first-principles calculations. They identified the migration barrier and the underlying mechanism. However, a complete study for the topology of the landscape around a stable minimum is still needed.

The principal goal of this study is to give a complete description of the diffusion mechanism in elementary semiconductors. As such, the SIEST-A-RT method is built up to develop an entire methodology for quantum mechanically precise and reasonably fast energy landscape exploration method. Our simulations are performed on supercells containing 215 atoms. This supercell size is computationally tractable and constitutes a reasonable size to simulate the host crystal. We generate the diffusion paths using the activation-relaxation technique (ART) Mousseau and Barkema (1998); Barkema and Mousseau (1996), which can sample efficiently the energy landscape of complex systems. The forces and energies are evaluated using SIESTA Sánchez-Portal et al. (1997); Soler et al. (2002), a self-consistent density functional method using standard norm-conserving pseudo-potentials and a flexible numerical linear combination of atomic orbital basis set. Combining these two methods allows us to identify very efficiently diffusion paths of various energy height.

This paper is constructed as follows. In the next section, we discuss the SIEST-A-RT method; simulation details are given in section III. We apply the method to the vacancy in silicon, the corresponding activated dynamics are detailed and discussed in section IV. Finally we emphasize the capability of the method in identifying competing events in silicon and its possible extention for future applications.

## II Method: SIEST-A-RT

In order to sample the diffusion paths associated with a vacancy in bulk Si, we integrate the activation-relaxation technique (ART nouveau) Barkema and Mousseau (1996); Malek and Mousseau (2000) to the ab-initio program SIESTA. This approach is similar to the integration of the hybrid eigenvector-following method to a tight-binding and an ab-initio code as performed by Munro, Kumeda and Wales Munro and Wales (1999); Kumeda et al. (2001) as well as the dimer-method Henkelman and Jónsson (1999) which was coupled with VASPKresse and Hafner (1993) by Henkelman and JònssonHenkelman et al. (2002).

As both ART nouveau and SIESTA are described elsewhere in the literature Malek and Mousseau (2000); Soler et al. (2002), we review these methods only briefly.

### II.1 ART nouveau

The activation-relaxation technique is a method for sampling efficiently the energy landscape of activated systems. Barkema and Mousseau (1996); Mousseau and Barkema (1998) Instead of following in details the thermal fluctuations, it searches directly for transition states; these are characterized as first-order saddle points in the energy landscape. This method was refined recently to ensure a controlled convergence to the saddle point (ART nouveau). Malek and Mousseau (2000) An ART event can be decomposed into the following steps. Activation: (1) Starting from a local minimum, push the configuration in a direction chosen randomly in the 3N-dimensional space; stop when the lowest eigenvalue becomes negative or when the configuration is far enough from the initial configuration. (2) Follow the eigendirection corresponding to this negative eigenvalue while minimizing the energy in the perpendicular directions; stop when the total force (parallel and perpendicular to the eigendirection) falls below a given threshold ($`0.1`$ eV/Å). This is the transition state. Relaxation: After pushing over the saddle point, minimize the energy using any standard minimization algorithm.

In order to sample only the activated mechanisms associated with the vacancy, we restrict the initial random projection for each event to the atoms surrounding the defect. Except for this initial selection, all atoms are allowed to move during the ART event.

ART uses the Lanczòs Lanczós (1988) algorithm to extract the lowest eigenvalue and its corresponding eigenvector; 16 to 20 force evaluations are sufficient, irrespective of the size of the system, to extract a stable eigenvalue value and eigenvector. Recent unpublished tests indicate that the number of force evaluations needed to reach a saddle point is about the same for the dimer method and ART nouveau Henkelman and Jónsson (1999), indicating no significant difference between the various activation-relaxation methods.

### II.2 Siesta

Forces and energies are evaluated using SIESTA Sánchez-Portal et al. (1997); Soler et al. (2002), a self-consistent density functional method using standard norm-conserving pseudo-potentials of Troullier-Martins Troullier and Martins (1991) factorized in the Kleiman-Baylander form Kleiman and Bylander (1982). Matrix elements are evaluated on a 3D grid in real space. The one-particle problem is solved using linear combination of pseudo-atomic orbitals (PAO) basis set of finite range. The quality of the basis can be improved by radial flexibilization, add more than one radial function within the same angular momentum (Multiple-$`\zeta `$), and angular flexibilization or polarization, add shells of different atomic symmetry (different l).

Because ART requires a numerical derivative of the force to compute the curvature of the energy landscape, it is essential to obtain highly converged forces. The various parameters used are discussed below.

### II.3 Optimization of the parameters

The introduction of a vacancy in the crystal affects considerably the relaxation of the simulation cell. Previous studies of the silicon vacancy using supercells ranging from 32 to 216 atoms with various Brillouin zone sampling techniques Puska et al. (1998); Kumeda et al. (2001); Virkkunen et al. (1993); Wang et al. (1991); Estreicher (2000); Puska (2000); Mercer et al. (1998); Probert and Payne (2003); Lento and Nieminen (2003) show a broad spread of formation energies (from 2.6 to 4.6 eV) and structural relaxation with ($`D_{2d}`$, $`C_{2v}`$, $`C_{3v}`$, $`T_d`$) symmetries. This scattering in mainly due to various convergence problems like basis set size, Brillouin zone sampling and supercell size. A recent work of Probert and Payne Probert and Payne (2003) propose a systematic methodology for accurate calculations of defect structure in supercells. We proceed in a similar way in using supercells with up to 512-atoms.

#### II.3.1 Mesh cut-off and the basis set

With the presence of a defect, the single-$`\zeta `$ basis (SZ), i.e., a single orbital per occupied state, is unsatisfactory: both the formation energy and reconstruction surrounding the vacancy are far from experimental value. Moreover, we find that the lack of overlap between the PAO induced discontinuities in the force, causing serious problems to the application of ART. It is possible to address these problems by using the optimized basis generated using the procedure of Anglada et al. Anglada et al. (2002) (SZ-optimized). Although much improved, the discontinuity in the force remains due to the presence of an underlying grid for the integration in real space which causes breaking of the translational symmetry. While the energy is converged with a real-space grid of about 0.03 Å (a mesh cut-off of 35-40 Ry), we find that we need a grid space of at most 0.02 Å (mesh cut-off of 50 Ry) to stabilize the force derivative; this spacing is comparable to the typical atomic displacement during the dynamics monitored by ART.

#### II.3.2 Brillouin zone sampling

Because the defect is highly localized even in large supercells, we need to have fully converged forces and total energies for the system. In order to isolate Brillouin zone sampling effects, we study force and energy convergence with respect to the density of k-points. We use a Monkhorst and Pack Monkhorst and Pack (1976) mesh to sample the Brillouin zone. Uniform meshes of $`q^3`$ shifted by ($`k_0,k_0,k_0`$) are used, with $`q`$ going from 1 to 4 and $`k_0`$ = 0 or 0.5. For $`q`$=1 no offsets are used ($`k_0=0`$) to sample the $`\mathrm{\Gamma }`$ point. Offsets of $`k_0=0.5`$ are used for $`q=2,3,4`$ for each cubic supercell size.

Figure 1 displays the total energy per atom embedded in differently sized supercells with respect to q. The total energy is converged to 1meV/atom at a mesh of $`3^3`$ for the 64 supercell, this corresponds to a density of 0.031 Å-1 calculated using our LDA lattice parameter 5.39 Å. As discussed by Probert and Payne Probert and Payne (2003) this density is sufficient to converge the total energy for larger systems. It is equivalent to a mesh of $`2^3`$ for 216 and 512 cells with densities of 0.031 Å-1 and 0.023 Å-1 respectively. Comparing the forces for unrelaxed vacancy using different meshes and supercells of 63, 215 and 511 atoms, we find that atomic forces converge faster than the total energy as a function of $`k`$-point sampling: convergence to 1 meV/Å per degree of freedom is reached with a $`2^3`$ mesh for the 63 cell and with the $`\mathrm{\Gamma }`$ point for larger cells.

#### II.3.3 Supercell calculations

As reported by Puska et al. Puska (2000), supercell method has obvious drawback because of the interaction of the defect and its periodic replicas. If the defect-defect distance is not large enough, the electronic structure of an isolated defect is distorted because the deep levels in the band gap form energy bands with a finite dispersion. The size of the supercell restricts also the ionic relaxation. The relaxation pattern is truncated midway between the defect and its nearest image.

In order to identify the smallest cell-size acceptable for the simulation of diffusion mechanisms with negligible size effects, we study the structural relaxation and the formation energies with respect to supercell size by employing supercells as large as 511 ionic sites.

Formation energies are calculated for fully converged supercells with respect to the basis set and Brillouin zone sampling as detailed in sections II.3.1II.3.2. Results for 63 and 215 atoms supercells are overestimated compared to the 511-cell results by 9 % and 2% respectively as shown in table 1 . This error on vacancy formation energy for small supercell was previously reported by Puska et al. Puska et al. (1998) as due to the interaction of the vacancy with its periodic images and reflects a repulsive component of the defect-defect interaction, this component tends to vanish for the 215 supercell.

Structural relaxation is also of a major importance, it informs us about the correctness of the atomic forces. The 63-atom supercell shows a wrong relaxation pattern. Tetrahedral ($`T_d`$ ) symmetry is a metastable minimum irrespective of the density of k points used with a formation energy of 3.99 eV ( see the second line on in table 1). Symmetry need to be broken artificially to get a Jahn-Teller distortion leading to a $`D_{2d}`$ symmetry with lower formation energy (3.59 eV). This is misleading and shows that size effects are important for supercells of 63 ionic sites and smaller. This size is insufficient to represent the system and the diffusion pattern correctly as discussed in more details in section IV.4.

For larger system, size effects tend to have minor effects on the relaxation geometry. Already at 215 cells, $`D_{2d}`$ symmetry is obtained even for low density of k-points ($`\mathrm{\Gamma }`$ only). The formation energy using $`\mathrm{\Gamma }`$ point for the 215 supercell is underestimated by 6% compared to the best converged results we got with the 511 supercell, it tends to converge to the formation energy of larger models relaxed with $`\mathrm{\Gamma }`$-point only.

We conclude, in agreement with previous works Probert and Payne (2003); Puska (2000); Puska et al. (1998) that a supercell of 215-atoms with $`\mathrm{\Gamma }`$-point is the minimal size and sampling we can use in order to ensure that both energy and forces are sufficiently accurate to provide results of a quality at least equal to experimental precision.

## III Details of the simulation

We perform first-principles electronic structure calculations based on the density functional theory within the local-density approximation (LDA) Ceperley and Alder (1980); Perdew and Zunger (1981) for the exchange-correlation energy. The static zero-temperature calculations for energies and forces are computed from SIESTA. Forces on the atoms are computed using the Hellman-Feynman theorem. The nuclear positions are optimized by using a conjugate gradient algorithm (CG) to minimize the total energy. Using the unit cell of two atoms and increasing the number of k points the equilibrium lattice constant of bulk Si converges to 5.39 Å. The cell volume and shape are kept fixed during the simulation. The crystalline 216 supercell is relaxed until the largest force component is of about 0.002 eV/Å. One conjugate gradient (CG) step is necessary to relax the structure with vanishingly small force components. To create the starting configuration with a defect, one atom is removed from the ideal crystalline cell. The obtained structure contains a mono-vacancy characterized by a vacancy-vacancy distance of 16.17 Å, with tetrahedral symmetry ($`T_d`$). The structure is then allowed to relax at constant volume. All the events are generated starting from a fully relaxed 215-atom unit-cell with a single vacancy displaying an expected Jahn-Teller distortion: the four atoms around the vacancy are paired two by two with $`D_{2d}`$ symmetry, as described in section IV.1. More than 120 events were started from this geometry, each launched in a different random direction, providing an extensive sampling of the energy landscape around the vacancy. The results from this search confirm a posteriori that this is indeed the minimum-energy configuration for the isolated vacancy. Using the set of parameters detailed in this section, one event with SIEST-A-RT requires about 700 force evaluations and takes about 3 days of CPU times on a single IBM Power 4 processor.

## IV Results

### IV.1 Neutral vacancy relaxation

The details of the relaxation around the vacancy depend on the local charge Watkins (1992); Bourgoin and Lanoo (1981). Since there is one electron per dangling bond for the neutral charge state ($`V^0`$), it is expected to undergo a Jahn-Teller distortion Watkins (1992); Bourgoin and Lanoo (1981) by forming two pairs out of the four dangling bonds of the vacancy. This is what we see. During the relaxation of the 215-atom cell, the total energy of the system is reduced by about 1.71 eV, due to the distortion of the lattice and the relaxation of the atoms around the vacancy. This relaxation energy is larger than the 0.36 eV as obtained from a study on 63-atom cell by Seong et Lewis Seong and Lewis (1996) and the 0.9 eV calculated on a HF cluster calculation Öğüt et al. (1997) and compared to values from a more recent study Probert and Payne (2003) 1.186 eV using GGA on a 256 supercell.

The relaxation can be described as follows. The vacancy nearest neighbors (NN) first move away from each other by pair along the $`<110>`$ axis Seong and Lewis (1996). The net local displacement of 0.39 Å per atom result into a tetragonal symmetry configuration ($`D_{2d}`$). The distorted tetrahedron formed by the atoms has four long bonds of 3.35 Å and two short bonds of 2.9 Å as shown in Fig. 3(A). The formation of bonds between the two pairs weakens the back bonds which are then significantly stretched Antonelli et al. (1998), up to 2.49 Å . The range of perturbation introduced by the defect in the 215 supercell affects up to 5 surrounding shells: the average amplitude of relaxation falls progessively from 0.39 Å for the first shell to 0.074 Å for the fifth shell which is of the order of atomic vibrations at room temperature.

The amplitude and the type of displacement compare very well with the previous results using DFT Puska et al. (1998); Kumeda et al. (2001); Virkkunen et al. (1993); Wang et al. (1991); Estreicher (2000); Puska (2000); Mercer et al. (1998); Probert and Payne (2003). The LDA results appear to be at odds with recent screened-exchange LDA calculations on partially relaxed neutral Si vacancy, it indicates that the relaxation could be outward Lento and Nieminen (2003). More work remains to establish if LDA really fails here. Our results on the pairing mode are compared with some of these results in table 1.

The formation energy $`E_f`$ is calculated by taking into account the distortion of the lattice, using the expression:

$$E_f=E_{N1}\frac{N1}{N}E_N$$
(1)
Where N is the number of atoms and $`E_N`$ is the relaxed total energy of the ideal structure. $`E_{N1}`$ is the total energy of the distorted system containing the defect. Within SIESTA, the formation energy of the mono-vacancy is 3.36 eV, which is in agreement with previous LDA and experimental studies. On the simulation side, Kelly et al. Kelly and Car (1992) calculated 3.5 eV, Blöchl et al. Blöchl et al. (1993) 4.1 eV, Mercer et al. Mercer et al. (1998) 3.63 eV and Kumeda et al. Kumeda et al. (2001) 3.32 eV. Experimentally, Watkins et al. Watkins and Corbett (1964) measured 3.6 $`\pm `$ 0.5 eV and Dannefaer et al. Dannefaer et al. (1986) 3.6$`\pm `$0.2 eV. A complete comparison in reported in table 1.

The energy of 511 supercell, $`E_f^{511}`$, is found to be 3.29 eV after a full convergence of the energy. Previous works by Wang et al. Wang et al. (1991) and Seong and Lewis Seong and Lewis (1996) reported $`E_f^{511}`$= 4.12 eV and 4.10, respectively. A recent study by Probert et al. give GGA results of a fully converged calculation for a vacancy in BCC silicon lattice with 256 atoms, they obtain a formation energy of 3.17 eV.

### IV.2 Simple diffusion paths

Even a defect as simple as a neutral vacancy can be associated with many activated mechanisms, some of which leading to migration: among the generated events, 60 % are simple diffusion mechanisms, while 15 % are reorientation and the remaining are more complex events involving a pair or more of active atoms.

#### IV.2.1 Nearest neighbor diffusion

The simplest migration process involves hopping of the vacancy to a nearest-neighbor site as shown in Fig. 2(a): a nearest-neighbor of the vacancy moves along the $`<111>`$ direction toward the vacant site (shown in Fig. 3A) and passes the metastable split-vacancy site before falling in the tetrahedral site previously formed by the vacancy. The energy barrier for this mechanism is $`0.40\pm 0.02`$ eV. The error estimate comes from the convergence criterion to the saddle point as explained in section II.1. This result for the migration energy is in good agreement with the experimental findings of Watkins et al. Watkins (1992) who measured $`0.45\pm `$0.04 eV. In a work similar to this one, but using the CPMD code with GGA (BLYP) functional, Kumeda et al. Kumeda et al. (2001)found a barrier of 0.58 eV.

The diffusion path is reconstructed in detail using the nudged-elastic band methodHenkelman et al. (2000) and relaxed until the total force is of 0.2 eV/Å. This path shows an unusually long transition plateau —ranging from $`d_r`$=0.7 to 1.2 Å— with a slight energy minimum at the split-vacancy site. At this point, the 3 back bonds are stretched and the diffusing atom forms weak bonds (length 2.685Å) with the surrounding atoms of the six membered ring. From this site, a small displacement of 0.11 Å is sufficient to overcome the barrier of 0.04 eV and complete the jump to the new minimum (see table 2). This value is of the order of the precision of the method; in addition, at room temperature this minimum would be smeared out by thermal vibrations. Because of a lack of a better approximation, we use Vineyard’s quasiharmonic approximations Vineyard (1957) to get a rough estimate of the attempt frequency (see section IV.2.3), approximating the transition plateau by a single barrier.

Although the total diffusion path in Fig. 2(a) is asymmetric, after separating the contribution coming from different moving atoms, we find that diffusion path is symmetric for the diffusing atom as shown in figure (Fig. 2(b)). The asymmetry in the path shown in Fig. 2(c) is due to the reconstruction: with four atoms reconstructing two-by-two, there are three possible orientations. Here, the orientation is different for the initial and final state. This should average out as the vacancy moves around.

A comparison of the formation energies using the Stillinger-Weber Maroudas and Brown (1993) and the MEAM Baskes et al. (1989) potentials for the vacancy and the split vacancy gives a migration energy of 0.43 and 0.37 eV, respectively. Although the formation energies for the vacancy are not well represented by empirical potentials the energy difference is in agreement with experiment and with our results. The minimum-energy path, however, is not well described by these potential and the local minima identified by Maroudas et al. Maroudas and Brown (1993) during the migration are absent in our case.

#### IV.2.2 Reorientation process

In addition to finding diffusion mechanisms associated with the vacancy, we also identify the trajectories responsible for atomic rearrangement in the reconstructed configuration. As there are three pairing possibilities for the $`D_{2d}`$ symmetry, it is possible for the configuration to hop from one of these states to either of the two others.

SIEST-A-RT generates easily this reorientation of the paired atoms. As illustrated in Fig. 4, starting with the reconstructed pairs formed initially — 1-4 and 2-3—, atom “3” moves away form its neighbor along the $`<111>`$ direction, breaks the weak bond with ”2” and reaches the transition state, which is in the trigonal symmetry ($`C_{3v}`$). The relaxation continues until new bonds form between the pairs 1-2 and 3-4. The activation energy for the reorientation of the neutral vacancy between two tetragonal orientations is 0.2 eV. Roberson et al. Roberson and Estreicher (1994) in a HF calculation scheme performed calculation on small clusters (44 silicon atoms) and found 0.37 eV; another cluster calculation, by Öğüt et al. Öğüt et al. (1997) found 0.32 eV.

The deformation to the $`C_{3v}`$ can be observed experimentally. The experiment consists of applying a stress at high temperature in the dark, then cooling down to about 20K, relieving the stress, illuminating the sample and monitoring the resulting vacancy species. The barrier measured from this alignment experiment is at 0.23 eV, in excellent agreement with our calculations Watkins (1992).

Our simulations also show that the pairing mode can change during the migration of the vacancy. For the simple diffusion plus reorientation, we find an overall energy barrier of 0.47 eV and a total displacement of 0.68 Å at the saddle point.

#### IV.2.3 Diffusion constants and diffusion rate

The diffusion is a two step process, consisting of the creation and the subsequent migration of a vacancy. Assuming that both processes are thermally activated, the activation energy $`E_a`$ for diffusion is the sum of the formation energy $`E_f`$ and the migration energy $`E_m`$. By collecting our data from the previous sections we get an activation energy of $`E_a`$ = 3.5 eV, which is in good agreement with $`E_a=4.14`$ eV recently measured by Bracht et al. Bracht et al. (1998), and with the activation enthalpy for self-diffusion $`H_V`$ = 4.07 $`\pm `$ 0.2 eV obtained by Tang et al. Tang et al. (1997) using tight-binding MD.

It is well known that there is a strong competition between the interstitial and the vacancy-mediated mechanisms of self-diffusion. These two contributions were separated in a metal-diffusion experiment Bracht et al. (1998) and the vacancy contribution to the self-diffusion was measured to be

$$C_V^{eq}(T)d_v(T)=0.92exp\frac{(4.14eV)}{k_BT}cm^2s^1$$
(2)
Other experimental data Tan and Gösele (1985) shows that

$$C_V^{eq}(T)d_v(T)=0.6exp\frac{(4.03eV)}{k_BT}cm^2s^1$$
(3)
Both data sets predict a self-diffusion dominated by vacancies at low temperature, where $`C_V^{eq}(T)`$ is the concentration of vacancies at a certain temperature.

Experiments under thermal equilibrium concentrations of native point defects have shown that the Si self-diffusion coefficient $`D_{Si}^{eq}`$ exhibits Arrhenius behavior,

$$D_{Si}^{eq}(T)=D_0e^{H_a/k_bT},$$
(4)
where $`H_a`$ is a single activation enthalpy, $`D_0`$ is the prefactor in a certain temperature range, $`T`$ denotes the absolute temperature, and $`k_B`$ is the Boltzman’s constant. Under the assumption that the diffusion proceeds through discrete jumps of equal length $`r`$ —the distance between nearest neighbor equilibrium sites, in diamond cubic lattices $`r=a\sqrt{3}/4`$, where $`a`$ is the lattice parameter — the prefactor is written as Agullo-Lopez et al. (1988)

$$D_0=\frac{z}{6}\mathrm{\Gamma }_0r^2$$
(5)
with $`\mathrm{\Gamma }_0`$ the jump, or attempt, frequency of the point defect from one equilibrium site to another, z denotes the number of neighbor sites, in the diamond lattice a vacancy can jump to one of the four neighbors. In the harmonic approximation, this attempt frequency is defined as Vineyard (1957); Phillips (2001)

$$\mathrm{\Gamma }_0=\frac{{\displaystyle \underset{i=1}{\overset{3N}{}}}\nu _i^{(i)}}{{\displaystyle \underset{i=1}{\overset{3N1}{}}}\nu _i^{(s)}}$$
(6)
where $`\nu _i^{(i)}`$ and $`\nu _i^{(S)}`$ are the normal mode frequencies at the minimum and the saddle point respectively and the product does not include the imaginary frequency at the saddle point. The eigenvalues are computed by diagonalizing the Hessian, obtained by numerical derivation with a step of 0.01 Å.

Accordingly, the migration entropy defined as the entropy difference between the saddle point and the equilibrium point configuration at low temperature can be written as

$$\mathrm{\Delta }S=S_V^m=k_B\mathrm{ln}\left[\frac{{\displaystyle \underset{i=1}{\overset{3N1}{}}}\nu _i^{(i)}}{{\displaystyle \underset{i=1}{\overset{3N1}{}}}\nu _i^{(s)}}\right]$$
(7)
The phonon frequency corresponding the direction of the jump is removed from the numerator. In addition the self diffusion entropy evaluated experimentally is the sum of the diffusion entropy and the formation entropy $`S_V^{SD}=S_V^f+S_V^m`$. Using the experimental $`S_V^{SD}=5.5k_B`$ by Bracht et al. Bracht et al. (1998) and a $`S_v^f=(5\pm 2)k_B`$ as calculated by Blöchl et al. Blöchl et al. (1993) and Dobson et al. Dobson et al. (1989), migration entropy is expected to lie within 1$`k_B`$ and 2.5$`k_B`$, our calculated data $`\mathrm{\Delta }S=3.036k_B`$.

The corresponding attempt frequency $`\mathrm{\Gamma }_0`$ = $`8.65\times 10^{12}s^1`$ is similar to the result of Maroudas et al. Maroudas and Brown (1993) ($`1.0539\times 10^{13}s^1`$), obtained with the Stillinger-Weber potential. The corresponding diffusion prefactor is $`D_0=3.141\times 10^3cm^2s^1`$. Our diffusion constant can be compared with results from Maroudas et al. Maroudas and Brown (1993) by taking into account the z factor that had been omitted, the resulting prefactor $`D_0=3.88\times 10^3cm^2s^1`$ is in agreement with our results. Tang et al. Tang et al. (1997); Jääskeläinen. et al. (2001) report $`D_0=1.18\times 10^4cm^2s^1`$ using TB-MD and thermodynamical integration, their migration energy ($`E_m`$) is of 0.1 eV.

Our numerical estimate of the diffusion rate $`\mathrm{\Gamma }=\mathrm{\Gamma }_0e^{E_m/k_BT}`$ at room temperature gives a characteristic time for diffusion ($`\tau =1/\mathrm{\Gamma }`$ ) of the order of microseconds. This timescale is not directly accessible to ab-initio molecular dynamics, so simulations must generally be run at temperatures close to the melting point. With SIEST-A-RT, we can study the activated mechanisms from the local energy minimum, allowing us to identify accurately subtle mechanisms such as that associated with the split-vacancy site.

### IV.3 Complex diffusion paths

The SIEST-A-RT sampling of the energy landscape around the vacancy also identifies a number of high-energy barrier events which involve a complex atomistic rearrangements. Two types of such events occur. The first one can be described as a diffusion event mediated by bond exchange. In this case, the overall result is a simple vacancy hop but with a more complex transition state and a barrier of 2.4 eV.

The second type of events involves what has been dubbed spectator events Kumeda et al. (2001) (see Fig. 5) where the diamond network surrounding the vacancy is alterated. These spectator events all correspond basically to the bond-switching mechanisms proposed many years ago by Wooten, Winer and Weaire as mechanism for amorphization Wooten et al. (1985). In recent years, the WWW mechanism has also been identified as the bond defect Cargnoni et al. (1998) and the four fold coordinated defect (FFCD) Goedecker et al. (2002).

Two types of spectator events are identified. In one of these events, two of the direct neighbors of the vacancy are pushed toward each other, leading to the saturation of the dangling bonds; as a consequence one of the long covalent bonds between the NN of the vacancy transform to a real bond in a form very similar to the bond defect complex identified by Marqués et al Marques et al. (2001). The lowest barrier we found is 2.63 eV for a process involving a direct neighbor of the vacancy (one of the bonds is already weak), the final configuration is 2.16 eV higher than the initial minimum. Kumeda et al. found the single spectator WWW barrier of 1.719 eV using LDA and 1.867 eV using BLYP functional Kumeda et al. (2001).

The second type, atoms far from the vacant site are involved. The effect of the vacancy on the WWW event barrier is less strong for complex events located far from the vacancy without involving a direct neighbor. As the lattice distortion increases, so does the energy barrier. For this mechanism, for example, the barrier ranges from 3.31-3.8 eV, a value which can be compared to the I-V recombination process identified to be responsible for the amorphization of crystalline silicon of Marqués et al Marques et al. (2001) the barrier is estimated to be between 4.58 and 4.75 eV using ab initio calculation Cargnoni et al. (1998). One can estimate the effect of the vacancy to lower the barrier by about 1.27 eV. Since the barrier is very high, this kind of events is unlikely to play any significant role in crystals at room temperature.

### IV.4 Size effects

As discussed, previously, it is necessary to use a large simulation cell in order to minimize the interaction between the images of the vacancy. This is particularly noticeable in a 63-atom supercell, even though the defect-defect separation is larger than 4 inter-atomic distances. It is this interaction that cause the unrealistic relaxation pattern identified in previous studies Mercer et al. (1998); Virkkunen et al. (1993), allowing the vacancy to relax both with a $`T_{2d}`$ and a $`D_{2d}`$ symmetries.

Relaxing a 63-atom supercell with SIESTA and the parameters described above, we find that two stable minima: the $`D_{2d}`$ (tetragonal) relaxation described previously as well as a $`T_d`$ (tetrahedral) configuration, a symmetry which is unstable in cells of 215 and 511 atoms. Within the later symmetry, all atoms relax inward, conserving the tetrahedral symmetry of the system at equal distance of 3.4 Å. The $`T_{2d}`$ structure can be reached with a single jump from the $`D_{2d}`$ site during a reorientation mechanism similar to that described in section IV.2.2, crossing a small barrier of of 94 meV with an overall displacement of 0.64 Å. The energy difference between the two symmetries is estimated to be 83 meV compared to 20 meV as obtained by Mercer et al. in a simulation also on a 63-atom cell Mercer et al. (1998). We obtain qualitatively similar results with $`2^3`$ k-points mesh.

The interaction between images also affects the barrier for the simple diffusion and the reorientation, which is underestimated by 50 % compared to the 215 cell results. Interestingly, this underestimation is compensated by a higher formation energy, leading to very similar activation energies for the two size systems.

The underestimation of the diffusion barrier in the 63-atom cell makes it easier for the vacancy to jump to the second-neighbor position. Although we saw many such events in the small cell, this mechanism was not generated in the 215-atom box; there does not appear to be any direct path for second-neighbor diffusion in silicon.

## V Discussion and Conclusions

In this paper, we report on the study of the activated dynamics associated with a neutral vacancy in bulk crystalline silicon. Following the work of Munro, Kumeda and Wales Munro and Wales (1999); Kumeda et al. (2001), we couple an ab initio code, SIESTA, with the activation-relaxation technique, ART nouveau, in order to sample the various activated mechanisms taking place around the vacancy accurately and efficiently. Defining moves directly in the energy landscape, it is as easy to generate the diffusion trajectory responsible for the high-energy bond-switching mechanisms than it is to find the path associated with the jump reorientation of the reconstructed state ( $`D_{2d}D_{2d}`$).

Simulating more than 120 activated trajectories, we find a number of different diffusion mechanisms associated with the vacancy. In particular, we recover the basic vacancy hopping diffusion mechanism and show that the path associated with it includes a very metastable state at the split-vacancy site. Moreover, we can reproduce without difficulty the jump between the symmetrically equivalent reconstructed states. We also generate a number of higher-energy diffusion mechanisms which reminiscent of those found in amorphous silicon as well as at the crystalline-amorphous interface.

Because we obtain, with good accuracy, the transition state associated with each of these mechanisms, we also compute the transition rate in the harmonic approximation. The result obtained is in good agreement with experiment as well as previous high-temperature TB-MD simulation Tang et al. (1997).

Overall, various barrier and relaxation energies computed here are in agreement with results obtained previously both experimentally and numerically. Using a large simulation cell and an extensive sampling of the energy landscape around the vacancy, we could also show that some of the previously obtained mechanisms were artifact due to the small size of the cell used. These results establish the validity of the SIEST-A-RT approach as a unique tool for the study of the diffusion and relaxation mechanisms of defects in crystalline materials. This method can be easily extended to study semiconductor compounds, adatom diffusion, interlayer diffusion, all systems with a landscape too complicated for high-temperature molecular dynamics. Even in system as simple as a vacancy in c-Si, we found a number of activated mechanism associated with non-trivial transition paths. The situation promises to be more complex in binary semiconductors, such as GaAs and GaN, where the bonding character and chemical composition are much more complex than for Si.

## VI Acknowledgments

This work is funded in part by NSERC and NATEQ. NM is a Cottrell Scholar of the Research Corporation. Most of the simulations were run on the computers of the Réseau québécois de calcul de haute performance (RQCHP) whose support is gratefully acknowledged.


text:
# Five nights of intensive R- and V-band photometry of QSO 0957+561A,B

## 1 Introduction

More than 20 years of monitoring of the ‘classical’ double quasar QSO 0957+561A,B (Walsh et al., 1979) has shown that this interesting object exhibits variability at several levels and time-scales. This variability has been extremely useful for such diverse applications as, e.g., estimating the Hubble constant via the ‘Refsdal time delay method’, constraining the continuum-emitting region of the quasar, and probing the mass distribution in the lens galaxy and its surrounding cluster. Despite its relative photometric simplicity, QSO 0957+561 ($`z_{\mathrm{QSO}}=1.41`$, $`z_{\mathrm{lens}}=0.36`$, $`\theta _{\mathrm{AB}}6.2`$ arcsec) has nevertheless proved more unwieldy than initially hoped. In this context we note that even the best sampled light curves may exhibit artificial fluctuations or seeming correlations, and the task of distinguishing between physical events (whether intrinsic or microlensing-induced) and spurious ones is not necessarily straightforward.

In this paper we employ our newly developed photometry scheme (Ovaldsen et al., 2003) and present two-colour photometry for QSO 0957+561A,B based upon about 1200 CCD images from five observing nights, four of which are consecutive. This densely sampled, high-quality data set allows us to, e.g., constrain the short time-scale variability and the noise level to be expected in past and future investigations on the same theme. Part of this study constitutes a re-reduction of existing data, and our results are at variance with previously published work (Colley et al., 2002, 2003, hereafter C02 and C03); for instance, we find no support for the hourly variability invoked by the above authors to facilitate a time delay determination with accuracy of a fraction of a day (417.09$`\pm `$0.07 d).

The data material is described in Sect. 2, while a short outline of the photometric procedure is given in Sect. 3. Photometry for comparison stars as well as for the A and B quasar components is presented and commented in Sect. 4. Sect. 5 concludes with a discussion and a summary.

## 2 Data set

Our original data set consists of 976 $`R`$-band and 229 $`V`$-band CCD frames obtained with the 2.56-m Nordic Optical Telescope (NOT) on La Palma, Canary Islands, Spain, using the ALFOSC instrument.111ALFOSC is owned by the Instituto de Astrofisica de Andalucia (IAA) and operated at the NOT under agreement between IAA and the NBIfAFG of the Astronomical Observatory of Copenhagen. The observations were made on 2000 January 25–26 and during the four nights 2001 March 14–18, and were part of the ‘QuOC Around the Clock’ 2000/2001 campaigns, designed to compare observing runs separated by the assumed time delay for QSO 0957+561 of $``$416 days; see C02 and C03, where a reduction of the $`R`$-band data from the same set is discussed. (Note that investigations in recent years seem to favour a lag around 425 days, see Ovaldsen et al., 2003 and references therein.) The detector is a 2k$`\times `$2k CCD with a gain of 1 e-/ADU and readout noise of approx. 6 e-/pixel. The nominal pixel size is 0.188 arcsec.

All images were pre-processed with bias subtraction and flat-fielding, exclusively by means of sky flats. The averaged flats typically represented a photon number 100–500 times the background level (which was between 100 and 600 photons) in the science frames. According to the online specification of the CCD chip, it contains 8 bad columns, 10 traps, and many dark ‘speckles’. Deplorably, we had to discard 107 frames ($`97R+10V`$) because one or both of the quasar components fell too close to a defective column.

The seeing varied from 0.6 to 2 arcsec (median value $`0.9`$). The point sources on the frames were usually slightly elliptical, with a median axis ratio of 1.1. However, some PSFs were highly elliptical, especially when the object was near the horizon.

## 3 Photometry scheme

The photometry package used to reduce the CCD images was developed by us (JT and JEO). The method is described in Ovaldsen et al. (2003), where it was applied to a five-year archival data set of QSO 0957+561; a full treatment is given in Ovaldsen (2002). Features include automatic source detection and localisation, field star aperture and PSF photometry, magnitude calibration, and aperture and PSF photometry of the twin quasar images (after lens galaxy subtraction and crosstalk correction).

### 3.1 Field star photometry and calibration

To obtain a reliable magnitude calibration, we believe it is imperative to include as many comparison stars as possible *and* to check their internal brightness consistency. We apply aperture photometry, using a 3 arcsec radius, to 7 comparison stars in the field (F, G, H, E, D, X, and R, see fig. 1 in Ovaldsen et al., 2003). The measured total intensities of the field stars are compared to the reference values (see below) to estimate the frame conversion factor. If more than one outlier is present, or if the conversion factor is statistically ambiguous, the frame is discarded. The obvious rationale for this policy is that if consistent photometry for the reference stars is impossible, it will be even more so for the quasar components.

Reference $`R`$-magnitudes for the comparison stars were kindly supplied by R. Schild (private communication). Lacking a proper set of reference $`V`$-magnitudes, we used the $`VR`$ colour terms of Serra-Ricart et al. (1999) for the stars H and D to convert our $`R`$-magnitudes for these two stars into $`V`$. Finally, the overall $`R`$\- and $`V`$-magnitude scales were adjusted to a best fit with the reference values.

It should be noted that the above procedure implies a large difference between our reference $`R`$-magnitudes and those of Serra-Ricart et al. We would like to stress, however, that this is not unusual. The reference magnitudes encountered in the literature (e.g. Keel, 1982, Vanderriest et al., 1989, and Serra-Ricart et al.), may differ by as much as 0.4 mag for a single star!

### 3.2 Photometry of the quasar images

The distinctive feature of QSO 0957+561 is the seeing-dependent contamination by the lensing galaxy, denoted G1, particularly of the closely juxtaposed B image. To correct for this influence, we construct and subtract from each frame a galaxy model. Using a composite $`R`$-band brightness profile derived by Bernstein et al. (1997) and the appropriate ellipse parameters, a seeing-free G1 model is synthesised and convolved with the PSF of the actual frame (see Ovaldsen, 2002 for details, especially as regards the magnitude calibration). The position and orientation is determined using the astrometry of Bernstein et al. and the calculated centroids of the quasar images. Note that the brightness level is adjusted by an amount estimated by a three-component (A, B, G1) simultaneous fitting in the entire data set. Without better options (see, however, Teuber et al., 2003), we use the same galaxy profile in the $`V`$-band, decreasing the brightness level accordingly.

After the galaxy subtraction, the PSF is used to correct for cross-contamination between the quasar components. Refraining from performing a PSF flux estimation as such, we have chosen the hybrid method consisting of a ‘PSF cleaning’ of one of the two quasar images followed by an aperture flux measurement of the other, and vice versa. If necessary, the procedure may be performed iteratively.

Separate corrections for galaxy contamination and aperture crosstalk are vital for achieving high-precision photometry, especially when there are large variations not only in seeing but also PSF ellipticity.

## 4 Photometric results

### 4.1 Field stars

Our calibration procedure accepted 817 (out of 879) $`R`$-band frames and 206 (out of 219) $`V`$-band frames. Photometry results for the comparison stars are summarised in Tables 1 and 2 for the $`R`$\- and $`V`$-band, respectively.

When presenting the quasar light curves (Sect. 4.3), we include the results for the fainter star, R, the brightness of which is comparable to that of the quasar components.

### 4.2 Galaxy contribution and aperture crosstalk

From the subtraction procedure described above, we found that, depending on seeing, the galaxy light in the $`R`$-band amounted to 18.24–18.26 and 20.45–20.25 mag in the B and A apertures, respectively. In the $`V`$-band, the corresponding values were 19.34–19.38 mag (B) and 21.50–21.30 mag (A). So, as seeing deteriorates, galaxy light spills out of the B aperture and into the A aperture. Crosstalk between the twin quasar images was less than 1 per cent for seeing values below 2 arcsec. Note that both corrections are functions also of PSF ellipticity, that is, axis ratio and position angle.

### 4.3 The A and B quasar light curves

Figs. 1 (2000 January), 2 and 3 (2001 March) display light curves of the A and B quasar images and the R star in the $`R`$\- and $`V`$-bands. No outliers have been removed. The nightly mean brightnesses of the same objects, along with the corresponding standard deviations, are given in Table 3.

Despite considerable hourly variations of such decisive observational parameters as airmass, seeing, PSF ellipticity, sky level, and S/N-ratio, the light curves exhibit no dramatic features and appear quite constant during each night. The nightly standard deviations for the quasar components range from 3.9 to 6.4 mmag, to be compared with the mean of the individual formal errors for each night of 4.6 to 7.1 mmag – in satisfactory agreement with the corresponding values for the reference R star. The formal uncertainties are rather conservative, as they include the errors from Poisson statistics, galaxy subtraction and magnitude calibration (see Ovaldsen, 2002 for details). We omit plotting error bars in the figures to increase clarity.

The most conspicuous features of the brightness data are the changes in A and B over a time span of the order of one year. It is found that A brightens by 0.12 mag in $`R`$ and 0.15 mag in $`V`$, whereas B becomes fainter by 0.043 mag in $`R`$ and 0.050 mag in $`V`$.

From Figs. 2 and 3 and Table 3, one might also deduce a small brightness decrease in A and B, in both colours, of 1–4 mmag/day. However, the individual data sets for each of the four consecutive nights reproduce little, if any, of this overall trend (see Sect. 5). For a discussion of the quasar variability over several weeks, see the recent paper by Ullán et al. (2003).

As the observing runs were separated by the approximate time delay of 416 days, we are now able to compare the brightness values of the A quasar component on 2000 January with the corresponding time-shifted ones of B (average of the 2001 March values). We find $`VR`$ colour indices of 0.357$`\pm `$0.007 mag for A and 0.318$`\pm `$0.009 mag for B. Converting the $`m_\mathrm{B}m_\mathrm{A}`$ differences from Table 3 into the physically significant B/A magnification ratios, we find values of 1.025$`\pm `$0.011 in $`R`$ and 1.062$`\pm `$0.009 in $`V`$. These results are in qualitative agreement with Serra-Ricart et al. (1999), but not quite consistent with the achromaticity of the gravitational lens and the expected differential reddening of B due to its larger optical depth (see also Sect. 5).

### 4.4 Comparison with previous reductions

As regards the quasar variability in the two periods, our results – typically constant light curves in both colours during each night – do not support the many short time-scale fluctuations seen in fig. 1 of C03 (and C02). We therefore suspect that most of these fluctuations may be spurious and a likely reflection of problems arising from combining photometry from 12 different observatories. Likewise, failure to reject images of poor quality – whether due to image flaws or, more subtly, to inconsistent reference star photometry – will inevitably introduce systematic errors which may easily be misidentified as brightness fluctuations.

The A and B brightness levels for both observing runs, 2000 January and 2001 March, also differ significantly from those of C02 and C03. As an example, we find that during the observing nights in the $`R`$-band 2001 March, A and B were on average 16.397 and 16.492 mag, respectively, i.e. A brighter than B. C03 (fig. 1) finds $`16.37`$ (A) and $`16.30`$ (B), i.e. B brighter than A. The above authors perform a single correction, solely a function of seeing, to account for both galaxy light contamination and aperture crosstalk. If, however, we *omit* the galaxy subtraction and only apply our crosstalk correction, we get average values for A and B in perfect agreement with C03. (Note that our corrections not only account for seeing effects on each particular frame, but also for PSF ellipticity.)

## 5 Discussion and summary

A careful reduction with our new photometry scheme has yielded quasar light curves with little noise and hardly any outliers among the more than 2000 A and B data points. The photometry program seems to cope well with the large nightly variations in airmass, seeing, PSF ellipticity, sky level, and S/N-ratio, and produces remarkably constant light curves for each night. Our light curves do not exhibit any hourly fluctuations, in contrast to what e.g. Colley & Schild (2000), C02 and C03 have found. As a consequence, it seems impossible to derive any firm time delay from the observational material discussed here, let alone to obtain an estimate accurate to five significant digits, as the authors of C03 have been able to.

The small decrease from one night to the next for both quasar images, in addition to the typically constant behaviour during each night, might be regarded as a tentative, short time-scale manifestation of the much-debated ‘zero lag correlation’ (between the theoretically unrelated A and B light curves) present in the brightness record of QSO 0957+561 (Kundic et al., 1995; Colley et al., 2003; Ovaldsen et al., 2003). In the last reference, several conspicuous, simultaneous fluctuations were found in the quasar light curves, with time-scales of up to several months and amplitudes up to 0.2 mag. However, four nights of continuous monitoring do admittedly not suffice to draw any firm conclusions about this issue. It constitutes nevertheless a strong indication that the correlation is a hitherto unexplained *observational* effect.

If not due to intrinsic fluctuations or microlensing, the observed small brightness variations – and the possible zero lag correlation as well – might be ascribed to colour effects arising from, e.g., changing observing conditions in connection with the spectral mismatch between quasar and comparison stars. However, the same remarks as in Sect. 4.3 regarding the large nightly changes in observational parameters coupled with the constant nightly light curves render this explanation unlikely. It remains a puzzle, though, that the night-to-night trend is not reflected in the nightly observations. Likewise, Ullán et al. (2003) found unexplained day-to-day variations that were very similar in both A and B. Because we have both constant and densely sampled light curves during the course of each night, we can exclude at least seeing, ellipticity, airmass and sky level from the list of possible explanations.

Combining our 2000 January A data with the 2001 March B data, in both colour bands, we have been able to estimate the optical continuum magnification ratio of QSO 0957+561A,B as well as the colour indices of the quasar images. The results stated in Sect. 4.3 seem somewhat surprising. There is no obvious explanation of the fact that A is redder than B; see however Serra-Ricart et al. (1999) and Michalitsianos et al. (1997). Next, the values of the flux ratio B/A differ in $`R`$ and $`V`$ (1.025 and 1.062, respectively). This might be an effect of inadequate galaxy modelling in the $`V`$-band (Sect. 3.2), but it seems very *ad hoc* to us to re-adjust the galaxy photometry for this reason only. For comparison in the $`R`$-band, we note that in the period 1992 October – 1995 June (true dates corresponding to the unshifted A image), the magnification ratio B/A was fairly stable at $``$1.05, based upon the $`R`$-band data set treated by Ovaldsen et al. (2003). The most straightforward explanation for this change is that one or both quasar components are affected by (long-term) microlensing.

Assuming that the observed gradient of, say, 3 mmag/day is typical for the unsystematic brightness variations in QSO 0957+561, whatever their origin, it follows from a simple random-walk argument that its expected variability over a few weeks will be of the order of 0.02 mag, in good agreement with the investigation by Ullán et al. (2003). Over a decade, the variations will be of the order of 0.2 mag, and such amplitudes are easily found in the brightness record (e.g. Pelt et al., 1998; Ovaldsen et al., 2003).

In conclusion, we hope that future observations and careful photometric analysis of the data from the elusive QSO 0957+561 system will shed light on some of the unresolved questions discussed above, not to mention the issues of time delay, lens modelling, and source structure.

## Acknowledgments

The work reported herein is based on observations made with the Nordic Optical Telescope, operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. We thank L. Goicoechea for exchanging ideas with us regarding variability and photometric problems relating to QSO 0957+561. The observations were made at the NOT by K. Kjernsmo, Ø. Saanum and A. K. D. Evans.


text:
# Jamming, two-fluid behavior and ‘self-filtration’ in concentrated particulate suspensions

## Abstract

We study the flow of model experimental hard sphere colloidal suspensions at high volume fraction $`\mathrm{\Phi }`$ driven through a constriction by a pressure gradient. Above a particle-size dependent limit $`\mathrm{\Phi }_0`$, direct microscopic observations demonstrate jamming and unjamming—conversion of fluid to solid and vice versa— during flow. We show that such a jamming flow produces a reduction in colloid concentration $`\mathrm{\Phi }_x`$ downstream of the constriction. We propose that this ‘self-filtration’ effect is the consequence of a combination of jamming of the *particulate* part of the system and continuing flow of the *liquid* part, i.e. the solvent, through the pores of the jammed solid. Thus we link the concept of jamming in colloidal and granular media with a ‘two-fluid’-like picture of the flow of concentrated suspensions. Results are also discussed in the light of Osborne Reynolds’ original experiments on dilation in granular materials.

In this paper we consider the pressure-driven flow of concentrated suspensions of model colloidal particles. Concentrated suspensions of particles in liquid solvents are ubiquitous in ‘soft matter’ technology (cosmetics, foods, building materials, paints, detergents, pharmaceuticals, waste management) as well as in natural phenomena (soil and wet sand, formation of porous rocks and sediments, landslip, etc). Much industrial processing of soft matter, and many examples of natural flow phenomena, involve pressure-driven flow; moreover the flow often features ‘complex’ geometries where convergent and divergent elements generate extensional components of strain (e.g. constricting and widening pipes in a transport system). Fundamental studies of suspension behavior in such complex flow remain rare in comparison to the wealth of studies of soft matter under simple shear Larson (1999). Rheometrical work on concentrated suspensions has demonstrated complicated effects such as stress-induced thickening, erratic flow response, and fluctuating viscosity Laun (2001); Frith et al. (1996).

In the rheology of very concentrated suspensions and other crowded soft matter systems a concept that has excited much recent speculation Liu and Nagel (1998) and theoretical and experimental work Bertrand et al. (2002); Holmes et al. (2003); Trappe et al. (2001); Cates et al. (1998); Liu and Nagel (2001) is that of *jamming*. Here we shall define jamming as *the conversion of a liquid system into a solid by imposed stress*. Jamming is very obvious in ‘hourglass’ flow of dry sand where stress-supporting solid arches form across the convergence, even though the sand typically flows more or less like a liquid in simpler geometries. ‘Dilation’ of wet sand is a related example with a venerable history, having been considered more than 100 years ago by Reynolds Reynolds (1885). However there remains no clear picture of the generic conditions required for nor the consequences of jamming in soft matter. Our study of a model system is aimed toward such a goal.

Our experimental system consists of polymethylmethacrylate (PMMA) spheres sterically stabilised by short grafted polymers (polyhydroxystearyl alcohol), suspended in a non-polar hydrocarbon solvent (decalin). This system is a good model of the conceptual ‘hard sphere’ suspension Pusey (1991). The particles can be manufactured in almost monodisperse batches (polydispersities typically less than 6-7%). Here we study three batches of particles of radii $`318\pm 10`$ nm, $`656\pm 20`$ nm and $`1000\pm 50`$ nm respectively, spanning the particle size range from ‘colloidal’ (where Brownian motion dominates in dilute suspensions) toward ‘granular’ where gravity is dominant.

We first demonstrate quantitatively one important and striking consequence of jamming in a concentrated suspension: that jamming in a convergent flow can induce a significant reduction in concentration downstream of the constriction, what one might call a ‘self-filtration’ behavior by the suspension. We study the effect of convergent flow on concentration by extracting a portion of a bulk sample through a constriction—put simply, by sucking a small volume out of a bulk cell through a narrow-barrelled syringe and comparing the resulting ‘extracted’ concentration $`\mathrm{\Phi }_x`$ with the bulk concentration $`\mathrm{\Phi }_b`$.

Samples are prepared as follows. Dilute suspensions are centrifuged until the particles form ‘solid’ sediments (see below for a discussion of jamming effects on sedimentation). To obtain given bulk volume fractions measured amounts of solvent are added to the sediments. In practice we prepare *pairs* of samples which have been subjected to exactly the same treatment (i.e. duration and rate of centrifugation etc. ) and which thus have equal $`\mathrm{\Phi }_b`$ to within the (small) uncertainties associated with measuring masses of sediment and added solvent. (We use *pairs* of samples for the comparison because, if $`\mathrm{\Phi }_x\mathrm{\Phi }_b`$, extracting some sample by syringe to measure $`\mathrm{\Phi }_x`$ and then measuring $`\mathrm{\Phi }_b`$ for the *same* sample would introduce systematic errors.) For one member of the pair we obtain the bulk colloidal *mass* fraction by scooping out a measured mass of the sample, allowing all solvent to evaporate, and measuring the mass of solids remaining. From the other member of the pair we extract a portion not by scooping but by *syringe*, and similarly measure the mass fraction of this extracted portion. We convert mass fractions into volume fractions $`\mathrm{\Phi }`$ (the more common measure of concentration in studies of colloidal suspensions) assuming a particle density $`\rho _p=1.188`$ g cm-3 (the bulk density of PMMA) and solvent density $`\rho _s=0.897`$ g cm-3. If there is no effect of convergent flow in the syringe on extracted volume fraction we should find $`\mathrm{\Phi }_x=\mathrm{\Phi }_b`$. A large number of experiments is carried out on rediluted and newly centrifuged pairs of samples to ensure that statistical variations due to measurement imprecision are minimised.

There are a number of ways of estimating $`\mathrm{\Phi }`$ for a colloidal suspension, none of which is without drawbacks. As is well known Pusey and van Megen (1986) the ‘mass fraction’ measure of $`\mathrm{\Phi }`$ tends to give systematically lower values compared to the other common measure in the PMMA system, that based on mapping experimental phase boundaries onto computer simulated fluid-crystal hard-sphere phase coexistence boundaries. Given the core-shell nature of the sterically-stabilised particles no simple measure can give an *absolutely* correct result in the sense of a ‘true hard sphere’ volume fraction. The important point for this work is that all measures are carried out in the same way so that we obtain volume fractions that are consistent and comparable with each other. In these experiments, given the small volumes (see below) of extracted suspension, the possible experimental errors are minimised by using the ‘mass fraction’ to obtain $`\mathrm{\Phi }`$; dilution and mapping onto phase boundaries is simply not possible with such small sample volumes. However care must be taken when comparing the numerical values of $`\mathrm{\Phi }`$ given here with those quoted elsewhere in the literature, which are often obtained by mapping onto the computer-simulation phase diagram.

All experiments are carried out in identical geometry: the ‘bulk’ cells are cylindrical glass cuvettes of diameter $`25`$ mm with sediment heights $`2.5`$ cm giving a total sample volume of $`12`$ cm3; the syringes used have maximum volume $`1.0`$ ml with an entry barrel internal diameter of $`1.6`$ mm; volumes extracted are typically $`0.40.5`$ ml, i.e. always $`<5\%`$ of the total bulk volume. The ‘sucking’ procedure is as follows: the syringe plunger is pulled out quickly by hand to generate an empty barrel volume of $`0.6`$ ml, so that a pressure drop ($`P_b/P_{atm}0.03`$) is instantaneously applied between the outside surface of the bulk at atmospheric pressure $`P_{atm}`$ and the surface inside the syringe barrel at reduced pressure $`P_b`$. In the majority of the experiments reported here we study this limiting case of a suddenly applied initial pressure drop (the pressure drop slowly decreases as the sample is extracted and enters the syringe barrel). A few experiments have also been carried out using larger initial pressure drop and using narrower constrictions.

Results for $`\mathrm{\Phi }_b`$ *vs* $`\mathrm{\Phi }_x`$ are shown in Figure 1, for the three particle sizes given above. As is clear, when $`\mathrm{\Phi }_b`$ is above some limit volume fraction $`\mathrm{\Phi }_0`$ the sample ’pumped’ through the contraction into the syringe barrel has a significantly reduced volume fraction $`\mathrm{\Phi }_x<\mathrm{\Phi }_b`$. The limit concentration $`\mathrm{\Phi }_0`$ is strongly dependent on particle size, decreasing with increasing particle size. Thus the effect is more visible with larger particles (approaching the granular scale) but even for the smallest particles, very much within the *colloidal* regime, there is a clear effect at the highest volume fractions. The reduction in extracted volume fraction is somewhat more severe with a smaller constriction geometry or with a larger initial pressure drop (see example open circles and open diamond, respectively, in Fig. 1), though we have yet to carry out detailed investigations on the effect of different pressure drop. In any case, since extraction of samples from bulk for e.g. observation in a microscope or measurement in a rheometer is very often achieved using syringes, pipettes etc., it is important for experimentalists to be aware of the sensitivity of key parameters such as volume fraction to such ‘processing’ prior to experiments. Of course similar flow situations involving concentrated suspensions are common in technological applications such as product delivery, waste processing, and so on.

Although it seems that this ‘self-filtration’ effect has not been reported in the colloidal literature, actually it is not too surprising if we consider the suspension as a ‘two-fluid’ system, that is a combination of *particulate fluid* and *liquid solvent*, where interaction between the particles may lead to a rheological ‘separation’ of the two fluids. There are a few studies of *squeeze flow* where a separation between solvent and dispersed phase is discussed Hebraud et al. (2000); Delhaye et al. (2000), but such a two-fluid picture of flow, though familiar in studies of polymer systems, is not often encountered in the colloidal literature. We propose that our results can be explained in a two-fluid spirit as an ‘extreme’ rheological separation of the particulate and solvent parts of the suspension. The particles *jam* and form a solid that resists the pressure drop and does not flow at all. Meanwhile the solvent, remaining liquid, cannot resist the pressure gradient and continues to flow through the pores of the jammed colloid. Hence there is an increased flow rate of solvent relative to particles, resulting in a reduced downstream particle volume fraction.

But is colloidal jamming really responsible for the measurements in Fig 1? To elucidate directly what is actually happening in our convergent pressure-driven flow we have carried out direct observations by optical microscope. To enable microscope observations, the PMMA particles are suspended in a mixture of solvents (decalin and tetralin) to partly match the particle refractive index and reduce multiple scattering, and samples are prepared (centrifuged) inside thin *rectangular* cuvettes. ‘Extraction’ is achieved by inserting a cylindrical glass capillary into the sediment, internal diameter $`1.0`$ mm, connected to a syringe whose plunger is withdrawn using a fixed-speed syringe pump. The flow in the region of the capillary tube entrance is observed with bright field microscopy, the field of view positioned at the region around the entrance to the capillary as shown in Fig. 2(d). Results in the form of digitised movies are available at jam , while figures 2(a)-(c) and 3 show extracted still-images. In all but one of the movies, the height of the field of view corresponds to $`1`$ mm and the entry to the capillary can be seen at the left of the picture. At such low magnification individual colloids cannot be resolved. In one movie (Ref. jam  experiment 4) the magnification is increased by a factor of 10 and the field of view is centred $`0.5`$ mm diagonally from the lower corner of the capillary entry; in this case the particles, $`1000`$ nm radius, are just resolvable.

Summarising our observations, flows at high $`\mathrm{\Phi }`$ are typically very erratic, demonstrating very clear *transient jamming* of the samples in the region of the tube entrance, involving slowing of flow, sudden ‘fracture’ events followed by speed-up, repeated re-jamming, and so on. Jamming as a transient ‘conversion’ of flowing liquid to stationary solid is strikingly apparent. In the movie at higher magnification (Ref. jam  experiment 4) the sample, blurred whilst flowing (i.e. liquid), can be seen to repeatedly momentarily freeze (become solid).

At the highest concentrations flow is localised near the entry of the capillary (Ref. jam  experiment 6). However dark lines form in arch-like shapes, often travelling in sudden ‘waves’ or ‘shocks’ out from the end of the capillary \[Figs. 2(b) and (c) and Fig 3\]. The ‘shocks’ are associated with subsequent speed-up of the flow, i.e. collapses of the jammed particle structure. Collapse is followed by rejamming, generating erratic changes in flow speed. The most likely optical origin of the dark lines is small-angle scattering (their appearance is very sensitive to refractive index contrast consistent with this), indicating localised changes in volume fraction. The exact nature of these waves or shocks of dark lines is unclear. They exhibit similarities with density or kinematic waves observed in granular media, such as ‘hourglass’ flow of sand Baxter et al. (1989). More detailed study and comparison with granular density waves will be the subject of future work Haw .

As $`\mathrm{\Phi }_b`$ is lowered toward the limit $`\mathrm{\Phi }_0`$, the jam-collapse-rejam behavior cycles faster (Ref. jam  experiments 3 and 4), suggesting that below some limiting concentration jams effectively collapse immediately, in other words the sample no longer jams during flow. Below $`\mathrm{\Phi }_0`$, we indeed observe smooth non-jamming flow of the suspension into the convergence \[Fig. 2(a) and Ref. jam  experiment 2\].

Continued flow of the reduced volume fraction (partially ‘self-filtered’) sample downstream presumably accounts for the *transience* of the jamming: jams give way just as, in hourglass flow of sand, solid arches give way to allow flow.

Collapse of jams apparently involves a sliding-solid fracture-like behavior (Ref. jam  experiment 5), sliding occurring at the dark lines visible in the sample. We speculate that such sliding-fracturing may also be a consequence of the combination of jamming and solvent permeation through the jam. Fractures may be localised by geometry-dependent pressure gradients driving solvent through the pores of the jams into regions of decreased particle concentration—sometimes called ‘microcracks’ in the rock fracture literature—which we associate with the dark lines we observe in the suspensions. The excess solvent in these regions lubricates localised sliding (‘cataclastic shear bands’ Main et al. (2000)) of opposing solid regions. Such a ‘lubricated slide’ picture has been proposed for slip of porous rocks in earthquakes (see e.g. Main et al. (2000)). The possibility that model suspensions might ultimately be used as simple ‘soft’ analogues of geological systems deserves further investigation.

A further demonstration of jamming is provided by simple observations of the behavior of sediments of the largest particles. Though apparently solid under careful handling, we have observed that the sediment formed after centrifugation very easily *reliquifies* under a slight deliberate *lateral* shaking. The sediment thus behaves as an example of what has been called ‘fragile matter’ Cates et al. (1998): it solidifies (jams) under application of a *unidirectional* force in the centrifuge (or over a longer time period in normal gravity) but cannot support (liquifies under) small stresses applied in *any other direction*.

Finally we note that the ‘sucking’ experiments may be compared to the early granular experiments of Osborne Reynolds, who squeezed rubber balloons filled with ball bearings and water Reynolds (1885). Reynolds observed that on applying pressure, water was drawn *into* the balloon—quite the reverse of the case of a balloon completely filled with a simple fluid, which would be forced *out*. The ball-bearing system inside the balloon can only strain by decreasing its local volume fraction, i.e. by *dilation* and increase of the total volume of the flexible balloon. Water is then naturally drawn in to fill the extra volume. Our observations in colloidal suspensions may also be interpreted in terms of dilation in the zone of convergent flow as follows. To allow flow in response to the applied pressure drop the region of sample just downstream of the convergence must dilate, i.e. decrease its volume fraction by taking on more solvent; but this solvent must come from somewhere, i.e. from the region upstream, hence generating an increase in upstream $`\mathrm{\Phi }`$. This increase in $`\mathrm{\Phi }`$ in turn leads to *jamming* of the colloidal particles.

The author acknowledges the support of the Royal Society of Edinburgh and thanks A. B. Schofield for the manufacture of the PMMA particles.