text:
# Density Estimators in Particle Hydrodynamics

## 1 Introduction

Smoothed Particle Hydrodynamics (SPH) has established itself as the workhorse for a variety of astrophysical fluid dynamical computations (Lucy Lc77 (1977), Ginghold & Monaghan GinMon77 (1977)). In a wide range of astrophysical environments this Lagrangian scheme offers substantial and often crucial advantages over Eulerian, usually grid-based, schemes. Astrophysical applications such as cosmic structure formation and galaxy formation, the dynamics of accretion disks and the formation of stars and planetary systems are testament to its versatility and succesful performance (for an enumeration of applications, and corresponding references, see e.g. the reviews by Monaghan Mon92 (1992) and Bertschinger Ber98 (1998)).

A crucial aspect of the SPH procedure concerns the proper estimation of the local density, i.e. the density at the location of the particles which are supposed to represent a fair – discrete – sampling of the underlying continuous density field. The basic feature of the SPH procedure for density estimation is based upon a convolution of the discrete particle distribution with a particular user-specified kernel function $`W`$. For a sample of $`N`$ particles, with masses $`m_j`$ and locations $`𝐫_j`$, the density $`\rho `$ at the location $`𝐫_i`$ of particle $`i`$ is given by

$$\rho (𝐫_i)=\underset{j=1}{\overset{N}{}}m_jW(𝐫_i𝐫_j,h_i),$$
(1)
in which the kernel resolution is determined through the smoothing scale $`h_i`$. Notice that generically the scale $`h_i`$ may be different for each individual particle, and thus may be set to adapt to the local particle density. Usually the functional dependence of the kernel $`W`$ is chosen to be spherically symmetric, so that it is a function of $`|𝐫_i𝐫_j|`$ only.

The evolution of the physical system under consideration is fully determined by the movement of the discrete particles. Given a properly defined density estimation procedure, the equations of motion for the set of particles are specified through a suitable Lagrangian, if necessary including additional viscous forces (see e.g. Rasio Ras99 (1999)).

In practical implementations, however, the SPH procedure involves a considerable number of artefacts. These stem from the fact that SPH particles represent functional averages over a certain Lagrangian volume. This averaging procedure is further aggravated by the fact that it is based upon a rather arbitrary user-specified choice of both the adopted resolution scale(s) $`h_i`$ and the functional form of the kernel $`W`$. Such a description of a physical system in terms of user-defined fuzzy clouds of matter is known to lead to considerable complications in realistic astrophysical circumstances. Often, these environments involve fluid flows exhibiting complex spatial patterns and geometries. In particular in configurations characterized by strong gradients in physical characteristics – of which the density, pressure and temperature discontinuities in and around shock waves represent the most frequently encountered example – SPH has been hindered by its relative inefficiency in resolving these gradients.

Given the necessity for the user to specify the characteristics and parameter values of the density estimation procedure, the accuracy and adaptibility of the resulting SPH implementation hinges on the ability to resolve steep density contrasts and the capacity to adapt itself to the geometry and morphology of the local matter distribution. A considerable improvement with respect to the early SPH implementations, which were based on a uniform smoothing length $`h`$, involves the use of adaptive smoothing lengths $`h_i`$ (Hernquist & Katz HK89 (1989)), which provides the SPH calculations with a larger dynamic range and higher spatial resolution. The mass distribution in many (astro)physical systems and circumstances is often characterized by the presence of salient anisotropic patterns, usually identified as filamentary or planar features. To deal with such configurations, additional modifications in a few sophisticated implementations attempted to replace the conventional – and often unrealistic and restrictive – spherically symmetric kernels by ones whose configuration is more akin to the shape of the local mass distribution. The corresponding results do indeed represent a strong argument for the importance of using geometrically adaptive density estimates. A noteworthy example is the introduction of ellipsoidal kernels by Shapiro et al. (SM96 (1996)). Their shapes are stretched in accordance with the local flow. Yet, while evidently being conceptually superior, their practical implementation does constitute a major obstacle and has prevented widescale use. This may be ascribed largely to the rapidly increasing number of degrees of freedom needed to specify and maintain the kernel properties during a simulation.

Even despite their obvious benefits and improvements, these methods are all dependent upon the artificial parametrization of the local spatial density distribution in terms of the smoothing kernels. Moreover, the specification of the information on the density distribution in terms of extra non physical variables, necessary for the definition and evolution of the properties of the smoothing kernels, is often cumbersome to implement and may introduce subtle errors (Hernquist H93 (1993), see however Nelson & Papaloizou NP94 (1994) and Springel & Hernquist SH02 (2002)). In many astrophysical applications this may lead to systematic artefacts in the outcome for the related physical phenomena. Within a cosmological context, for example, the X-ray visibility of clusters of galaxies is sensitively dependent upon the value of the local density, setting the intensity of the emitted X-ray emission by the hot intergalactic gas. This will be even more critical in the presence of feedback processes, which for sure will be playing a role when addressing the amount of predicted star formation in simulation studies of galaxy formation.

Here, we seek to circumvent the complications induced by the kernel parametrization and introduce and propose an alternative to the use of kernels for the quantification of the density within the SPH formalism. This new method, based upon the Delaunay Tessellation Field Estimator (DTFE, Schaap & van de Weygaert SW00 (2000)), has been devised to mould and fully adapt *itself* to the configuration of the particle distribution. Unlike conventional SPH methods, it is able to deal self-consistently and naturally with anisotropies in the matter distribution, even when it concerns caustic-like transitions. In addition, it manages to succesfully treat density fields marked by structural features over a vast (dynamic) range of scales.

The DTFE produces density estimates on the basis of the particle distribution, which is supposed to form a discrete spatial sampling of the underlying continuous density field. As a linear multidimensional field interpolation algorithm it may be regarded as a first-order version of the natural neighbour algorithm for spatial interpolation (Sibson Sib81 (1981), also see e.g. Okabe et al. Oketal00 (2000)). In general, applications of the DTFE to spatial point distributions have demonstrated its success in dealing with the complications of anisotropic geometry and dynamic range (Schaap & van de Weygaert SW00 (2000)). The key ingredient of the DTFE procedure is that of the Delaunay triangulation, serving as the complete covering of a sample volume by mutually exclusive multidimensional linear interpolation intervals.

Delaunay tessellations (Delone 1934; see e.g. Okabe et al. 2000 for extensive review) form the natural framework in which to discuss the properties of discrete point sets, and thus also of discrete samplings of continuous fields. Their versatility and significance have been underlined by their widespread applications in such areas as computer graphics, geographical mapping and medical imaging. Also, they have already found widespread application in a variety of ‘conventional’grid-based fluid dynamical computation schemes. This may concern their use as a non-regular application-oriented grid covering of physicalsystems, which represents a prominent procedure in technological applications. More innovating has been their use in Lagrangian ‘moving-grid’ implementations (see Mavripilis Mv97 (1997) for a review, and Whitehurst Wh95 (1995) for a promising astrophysical application).

It seems therefore a good idea to explore the possibilities of applying the DTFE in the context of a numerical hydrodynamics code. Here, as a first step, we wish to obtain an idea of the performance of a hydro code involving the use of DTFE estimates with respect to an equivalent code involving regular SPH density estimates. The quality of the new DTFE method with respect to the conventional SPH estimates, and their advantages and disadvantages under various circumstances, are evaluated by a comparison between the density field which would be yielded by a DTFE processing of the resulting SPH particle distribution and that of the regular SPH procedure itself. In this study, we operate along these lines by a comparison of the resulting matter distributions in the situation of a representative stochastic multiphase density field. This allows us to make a comparison between both density estimates in a regime for which an improved method for density estimates would be of great value. We should point out a major drawback of our approach, in that we do not really treat the DTFE density estimate in a self-consistent fashion. Instead of being part of the dynamical equations themselves we only use it as an analysis tool of the produced particle distribution. Nevertheless, it will still show the value of the DTFE in particle gasdynamics and give an indication of what kind of differences may be expected when incorporating in a fully self-consistent manner the DTFE estimate in an hydrocode.

On the basis of our study, we will elaborate on the potential benefits of a hydrodynamics scheme based on the DTFE. Specifically, we outline how we would set out to develop a complete particle hydrodynamics code whose artificial kernel based nature is replaced by the more natural and self-adaptive approach of the DTFE. Such a DTFE based particle hydrodynamics code would form a promising step towards the development of a fully tessellation based quasi-Eulerian moving-grid hydrodynamical code. Such would yield a major and significant step towards defining a much needed alternative and complement to currently available simulation tools.

## 2 DTFE and SPH density estimates

The methods we use for SPH and DTFE density estimates have been extensively described elsewhere (Hernquist & Katz HK89 (1989), Schaap & van de Weygaert SW00 (2000)). Here, we will only summarize their main, and relevant, aspects.

### 2.1 SPH density estimate

Amongst the various density recepies employed within available SPH codes, we use the Hernquist & Katz (HK89 (1989)) symmetrized form of Eq. 1, using adaptive smoothing lengths:

$$\widehat{\rho }_i=\frac{1}{2}\underset{j}{}m_j\left\{W(|r_ir_j|,h_i)+W(|r_ir_j|,h_j)\right\},$$
(2)
The smoothing lengths $`h_i`$ are chosen such that the sum involves around 40 nearest neighbours. For the kernel $`W`$ we take the conventional spline kernel described by Monaghan (Mon92 (1992)). Other variants of the SPH estimate produce comparable results.

### 2.2 DTFE density estimate

The DTFE density estimating procedure consists of three basic steps.

Starting from the sample of particle locations, the first step involves the computation of the corresponding Delaunay tessellation. Each Delaunay cell $`T_m`$ is the uniquely defined tetrahedron whose four vertices (in 3-D) are the set of 4 sample particles whose circumscribing sphere does not contain any of the other particles in the set. The Delaunay tessellation is the full covering of space by the complete set of these mutually disjunct tetrahedra. Delaunay tessellations are well known concepts in stochastic and computational geometry (Delaunay Del34 (1934), for further references see e.g. Okabe et al. Oketal00 (2000), Møller Mol94 (1994) and van de Weygaert Wey91 (1991)).

The second step involves estimating the density at the location of each of the particles in the sample. From the definition of the Delaunay tessellation, it may be evident that there is a close relationship between the volume of a Delaunay tetrahedron and the local density of the generating point process (telling examples of this may be seen in e.g. Schaap & van de Weygaert 2002a). Evidently, the “empty” cirumscribing spheres corresponding to the Delaunay tetrahedra, and the volumes of the resulting Delaunay tetrahedra, will be smaller as the number density of sample points increases, and vice versa. Following this observation, a proper density estimate $`\widehat{\rho }`$ at the location $`𝐱_i`$ of a sampling point $`i`$ is obtained by determining the properly calibrated inverse of the volume $`𝒲_{\mathrm{Vor},i}`$ of the corresponding contiguous Voronoi cell. The contiguous Voronoi cell $`𝒲_{\mathrm{Vor},i}`$ is the union of all Delaunay tetrahedra $`T_{m,i}`$ of which the particle $`i`$ forms one of the four vertices, i.e. $`𝒲_{\mathrm{Vor},i}=_mT_{m,i}`$. In general, when a particle $`i`$ is surrounded by $`N_T`$ Delaunay tetrahedra, each with a volume $`𝒱(T_{m,i})`$, the volume of the resulting contiguous Voronoi cell is

$$𝒲_{\mathrm{Vor},i}=\underset{m=1}{\overset{N_T}{}}𝒱(T_{m,i}).$$
(3)
Note that $`N_T`$ is not a constant, but in general may acquire a different value for each point in the sample. For a Poisson distribution of particles this is a non-integer number in the order of $`N_T27`$ (van de Weygaert Wey94 (1994)). Generalizing to an arbitrary $`D`$-dimensional space, and assuming that each particle $`i`$ has been assigned a mass $`m_i`$, the estimated density $`\widehat{\rho }_i`$ at the location of particle $`i`$ is given by (see Schaap & van de Weygaert SW00 (2000))

$$\widehat{\rho }(𝐫_i)=(D+1)\frac{m_i}{𝒲_{\mathrm{Vor},i}},$$
(4)
In this, we explicitly express $`𝒲_{\mathrm{Vor},i}`$ for the general $`D`$-dimensional case. The factor $`(D+1)`$ is a normalization factor, accounting for the $`(D+1)`$ different contiguous Voronoi hypercells to which each Delaunay hyper“tetrahedron” is assigned, one for each vertex of a Delaunay hyper“tetrahedron”.

The third step is the interpolation of the estimated densities $`\widehat{\rho }_i`$ over the full sample volume. In this, the DTFE bases itself upon the fact that each Delaunay tetrahedron may be considered the natural multidimensional equivalent of a linear interpolation interval (see e.g Bernardeau & van de Weygaert BerWey96 (1996)). Given the $`(D+1)`$ vertices of a Delaunay tetrahedron with corresponding density estimates $`\widehat{\rho }_j`$, the value $`\widehat{\rho }(𝐫)`$ at any location $`𝐫`$ within the tetrahedron can be straightforwardly determined by simple linear interpolation,

$$\widehat{\rho }(𝐫)=\widehat{\rho }(𝐫_{i0})+(\widehat{\rho })_{\mathrm{Del},m}(𝐫𝐫_{i0}),$$
(5)
in which $`𝐫_{i0}`$ is the location of one of the Delaunay vertices $`i`$. This is a trivial evaluation once the value of the (linear) density gradient $`(\widehat{\rho })_{\mathrm{Del},m}`$ has been estimated. For each Delaunay tetrahedron $`T_m`$ this is accomplished by solving the the system of $`D`$ linear equations corresponding to each of the remaining $`D`$ Delaunay vertices constituting the Delaunay tetrahedron $`T_m`$. The “minimum triangulation” property of Delaunay tessellations underlying this linear interpolation, minimum in the sense of representing a volume-covering network of optimally compact multidimensional “triangles”, has been a well-known property utilized in a variety of imaging and surface rendering applications such as geographical mapping and various computer imaging algorithms.

### 2.3 Comparison

Comparing the two methods, we see that in the case of SPH the particle ‘size’ and ‘shape’ (i.e. its domain of influence) is determined by some arbitrary kernel $`W(r,h_i)`$ and a fortuitous choice of smoothing length $`h_i`$ (assuming, along with the major share of SPH procedures, a radially symmetric kernel). In the case of the DTFE method the particles’ influence region is fully determined by the sizes and shapes of the Delaunay cells $`T_{m,i}`$, themselves solely dependent on the particle distribution. In other words, in regular SPH the density is determined through the kernel function $`W(𝐱)`$, while in DTFE it is solely the particle distribution itself setting the estimated values of the density. Contrary to the generic situation for the kernel dependent methods, there are no extra variables left to be determined. One major additional advantage is that it is therefore not necessary to worry about the evolution of the kernel parameters.

Both methods do display some characteristic artefacts in their density reconstructions (see Fig. 1). To a large extent these may be traced back to the implicit assumptions involved in the interpolation procedures, a necessary consequence of the finite amount of information contained in a discrete representation of a continuous field. SPH density fields implicitly contain the imprint of the specified and applied kernel which, as has been discussed before, may seriously impart its resolving power and capacity to trace the true geometry of structures. The DTFE technique, on the other hand, does produce triangular artefacts. At instances conspicuously visible in the DTFE reconstructed density fields, they are the result of the linear interpolation scheme employed for the density estimation at the locations not coinciding with the particle positions. In principle, this may be substantially improved by the use of higher order interpolation schemes. Such higher-order schemes have indeed been developed, and the ones based upon the natural neighbour interpolation prescription of Sibson (Sib81 (1981)) have already been succesfully applied to two-dimensional problems in the field of geophysics (Sambridge et al. Sametal95 (1995), Braun & Sambridge BraSam95 (1995)) and solid state physics (Sukumar Suk98 (1998)).

## 3 Case study: two-phase interstellar medium

For the sake of testing and comparing the SPH and DTFE methods, we assess a snapshot from a simulation of the neutral ISM. The model of the ISM is chosen as an illustration rather than as a realistic model.

The “simulation” sample of the ISM consists of HI gas confined in a periodic simulation box with a size $`L=0.6`$kpc3. The initially uniform density of the gas is $`n_H=0.3`$ cm-3, while its temperature is taken to be $`T=10000`$K. No fluctuation spectrum is imposed to set the initial featureless spatial gas distribution. To set the corresponding initial spatial distribution of the $`N=64000`$ simulation particles, we start from relaxed initial conditions according to a “glass” distribution (e.g. White white94 (1994)).

The evolution of the gas is solely a consequence of fluid dynamical and thermodynamical processes. No self gravity is included. As for the thermodynamical state of the gas, cooling is implemented using a fit to the Dalgarno-McCray (DalCra72 (1972)) cooling curve. The heating of the gas is accomplished through photo-electric grain heating, attributed to a constant FUV background ($`1.7`$ G0, with G0 the Habing field) radiation field. The parameters are chosen such that after about 15 Myrs a two-phase medium forms which consists of warm ($`10000`$ K) and cold ($`>100`$ K) HI gas.

The stage at which a two-phase medium emerges forms a suitable point to investigate the performance of the SPH and DTFE methods. At this stage we took a snapshot from the simulation, and subjected it to further analysis. For a variety of reasons, the spatial gas distribution of the snapshot is expected to represent a challenging configuration. The multiphase character of the resulting particle configuration is likely to present a problem for regular SPH. Density contrasts of about four orders of magnitude separate dense clumps from the surrounding diffuse medium through which they are dispersed. Note that a failure to recover the correct density may have serious repercussions for the computed effects of cooling. In addition, we notice the presence of physical structures with conspicuous, aspherical geometries (see Fig. 1 & 2), such as anisotropic sheets and filaments as well as dense and compact clumps, which certainly do form a challenging aspect for the different methods.

### 3.1 Results

Fig. 1 offers a visual impression of the differences in performance between the SPH and DTFE density reconstructions. The greyscale density maps in Fig. 1 (lower left: SPH, lower right: DTFE) represent 2-D cuts through the corresponding 3-D density field reconstructions (note that contrary to the finite width of the corresponding particle slice, upper left frame, these constitute planes with zero thickness).

Immediately visible is the more crispy appearance of the DTFE density field, displaying substantially more contrast in conjunction with more pronounced structural features. Look e.g. at the compact clump in the lower righthand corner ($`X0.5,Y0.12`$), forming a prominent and tight spot in the DTFE density field. The clump at ($`X0.48,Y0.52`$) represents another telling example, visible as a striking peak in the DTFE rendering while hardly noticeable in the SPH reconstruction. Structures in the SPH field have a more extended appearance than their counterparts in the DTFE field, whose matter content has been smeared out more evenly, over a larger volume, yielding features with a significantly lower contrast. In this assessment it becomes clear that the DTFE reconstruction adheres considerably closer to the original particle distribution (top lefthand frame). Apparently the DTFE succeeds better in rendering the shapes, the coherence and the internal composition in the displayed particle distribution. At various locations, the DTFE even manages to capture structural details which seem to be absent in the SPH density field.

To quantify the visual impressions of Fig. 1, and to analyze the nature of the differences between the two methods, we plot the ratio $`\rho _{\mathrm{DTFE}}/\rho _{\mathrm{SPH}}`$ as a function of the SPH density estimate $`\rho _{\mathrm{SPH}}/\rho `$ (in units of the average density $`\rho `$). Doing so for all particles in the sample (Fig. 1, top righthand, Fig. 2, top lefthand) immediately reveals interesting behaviour. The scatter diagram does show that the discrepancies between the two methods may be substantial, with density estimates at various instances differing by a factor of 5 or more.

Most interesting is the finding that we may distinguish clearly identifiable and distinct regimes in the scatter diagram of $`\rho _{\mathrm{DTFE}}/\rho _{\mathrm{SPH}}`$ versus $`\rho _{\mathrm{SPH}}/\rho `$. Four different sectors may be identified in the scatter diagram. Allowing for some arbitrariness in their definition, and indicating these regions by digits 1 to 4, we may organize the particles according to density-related criteria, roughly specified as (we refer to Fig. 2, top left frame, for the precise definitions of the domains):

1. low density regions:
$`\rho _{\mathrm{SPH}}/\rho <\mathrm{\hspace{0.17em}1}`$
2. medium density regions, DTFE smaller than SPH:
$`\rho _{\mathrm{DTFE}}<\rho _{\mathrm{SPH}}`$; $`1<\rho _{\mathrm{SPH}}/\rho <10`$
3. medium density regions, DTFE larger than SPH:
$`\rho _{\mathrm{DTFE}}>\rho _{\mathrm{SPH}}`$; $`1<\rho _{\mathrm{SPH}}/\rho <10`$
4. high density regions:
$`\rho _{\mathrm{DTFE}}\rho _{\mathrm{SPH}}`$; $`\rho _{\mathrm{SPH}}/\rho >\mathrm{\hspace{0.17em}10}`$

The physical meaning of the distinct sectors in the scatter diagram becomes apparent when relating the various regimes with the spatial distribution of the corresponding particles. This may be appreciated from the five subsequent frames in Fig. 2, each depicting the related particle distribution in the same slice of width 0.04 kpc. The centre and bottom frames, numbered 1 to 4, show the spatial distribution of each group of particles, isolated from the complete distribution (top right frame, Fig. 2). These particle slices immediately reveal the close correspondence between any of the sectors in the scatter diagram and typical features in the spatial matter distribution of the two-phase interstellar medium. This systematic behaviour seems to point to truly fundamental differences in the workings of the SPH and DTFE methods, and would be hard to understand in terms of random errors. The separate spatial features in the gas distribution seem to react differently to the use of the DTFE method.

We argue that the major share of the disparity between the SPH and DTFE density estimates has to be attributed to SPH, mainly on the grounds of the known fact that SPH is poor in handling nontrivial configurations such as encountered in multiphase media. By separately assessing each regime, we may come to appreciate how these differences arise. In sector 1, involving the diffuse low density medium, the DTFE and SPH estimates are of comparable magnitude, be it that we do observe a systematic tendency. In the lowest density realms, whose relatively smooth density does not raise serious obstacles for either method, DTFE and SPH are indeed equal (with the exception of variations to be attributed to random noise). However, near the edges of the low density regions, SPH starts to overestimate the local density as the kernels do include particles within the surrounding high density structures. The geometric interpolation of the DTFE manages to avoid this systematic effect (see e.g. Schaap & van de Weygaert 2002a,b), which explains the systematic linear decrease of the ratio $`\rho _{\mathrm{DTFE}}/\rho _{\mathrm{SPH}}`$ with increasing $`\rho _{\mathrm{SPH}}/\rho `$. To the other extreme, the high density regions in sector 4 are identified with compact dense clumps as well as with their extensions into connecting filaments and walls. On average DTFE yields higher density estimates than SPH, frequently displaying superior spatial resolution (see also greyscale plot in Fig. 1). Note that the repercussions may be far-reaching in the context of a wide variety of astrophysical environments characterized by strongly density dependent physical phenomena and processes ! The intermediate regime of sectors 2 and 3 clearly connects to the filamentary structures in the gas distribution. Sector 2, in which the DTFE estimates are larger than those of SPH, appears to select out the inner parts of the filaments and walls. By contrast, the higher values for the SPH produced densities in sector 3 are related to the outer realms of these features. This characteristic distinction can be traced back to the failure of the SPH procedure to cope with highly anisotropic particle configurations. While it attempts to maintain a fixed number of neighbours within a spherical kernel, it smears out the density in a direction perpendicular to the filament. This produces lower estimates in the central parts, which are compensated for with higher estimates in the periphery. Evidently, the adaptive nature of DTFE does not appear to produce similar deficiencies.

## 4 The DTFE particle method

Having demonstrated the improvement in quality of the DTFE density estimates, this suggests a considerable potential for incorporating the DTFE in a self-consistent manner within a hydrodynamical code. Here, we first wish to indicate a possible route for accomplishing this in a particle hydrodynamics code through replacement of the kernel based density estimates (1) by the DTFE density estimates. We are currently in the process of implementing this. The formalism on which this implementation is based can be easily derived, involving nontrivial yet minor modifications. Essentially, it uses the same dynamic equations for gas particles as those in the regular SPH formalism, the fundamental adjustment being the insertion of the DTFE densities instead of the regular SPH ones. In addition, a further difference may be introduced through a change in treatment of viscous forces. Ultimately, this will work out into different equations of motion for the gas particles. A fundamental property of a DTFE based hydrocode, by construction, is that it conserves mass exactly and therefore obeys the continuity equation. This is not necessarily true for SPH implementations (Hernquist & Katz HK89 (1989)).

The start of the suggested DTFE particle method is formed by the discretized expression for the Lagrangian $`L`$ for a compressible, nondissipative flow,

$$L=\underset{i}{}m_i\left(\frac{1}{2}v_i^2+u_i(\rho _i,s_i)\right),$$
(6)
where $`m_i`$ is the mass of particle $`i`$, $`v_i`$ its velocity, $`s_i`$ the corresponding entropy and $`u_i`$ its specific internal energy. In this expression, $`\rho _i`$ is the density at location $`i`$, as yet unspecified. The resulting Euler-Lagrange equations are

$$\frac{dv_i}{dt}=\underset{j}{}m_j\left(\frac{u_j}{\rho _j}\right)_s\frac{\rho _j}{x_i}.$$
(7)
The standard SPH equations of motion then follow after inserting the SPH density estimate (Eq. 1). Instead, insertion of the DTFE density (Eq. 4) will lead to the corresponding equations of motion for the DTFE-based formalism. Note that the usual conservation properties related to Eq. 6 remain intact. After some algebraic manipulation, thereby using the basic thermodynamic relation for a gas with equation of state $`P(\rho )`$,

$$\left(\frac{u_i}{\rho _i}\right)_s=\frac{P_i}{\rho _i^2},$$
(8)
we finally obtain the equations of motion for the gas particles (moving in $`D`$-dimensional space),

$$\frac{dv_i}{dt}=\frac{1}{D+1}\underset{m=1}{\overset{N_T}{}}P(T_{m,i})\frac{𝒱(T_{m,i})}{x_i}.$$
(9)
This expression involves a summation over all $`N_T`$ Delaunay tetrahedra $`T_{m,i}`$, with volumes $`𝒱(T_{m,i})`$, which have the particle $`i`$ as one of its four vertices. The pressure term $`P(T_{m,i})`$ is the sum over the pressures $`P_j`$ at the four vertices $`j`$ of tetrahedron $`T_{m,i}`$, $`P(T_{m,i})=P_j`$.

As an interesting aside, we point out that unlike in the conventional SPH formalism, this procedure implies an exactly vanishing acceleration $`dv_i/dt`$ in the case of a constant pressure $`P`$ at each of the vertices of the Delaunay tetrahedra containing particle $`i`$ as one of their vertices. The reason for this is that one can then invoke the definition of the volume of the contiguous Voronoi cell corresponding to point $`i`$ (Eqn. 3), yielding

$$\frac{dv_i}{dt}=\frac{1}{D+1}P\frac{𝒲_{Vor,i}}{x_i}.$$
(10)
Since the volume of the contiguous Voronoi cell does not depend on the position of particle $`i`$ itself (it lies in the interior of the contiguous Voronoi cell), the resulting acceleration vanishes. Another interesting notion, which was pointed out by Icke (ick2002 (2002)), is that Delaunay tessellations also provide a unique opportunity to include a natural treatment of the viscous stresses in the physical system. We intend to elaborate on this possibility in subsequent work dealing with the practical implementation along the lines sketched above.

## 5 Delaunay tessellations
and ‘moving grid’ hydrocodes

Ultimately, the ideal hydrodynamical code would combine the advantages of the Eulerian as well as of the Lagrangian approach. In their simplest formulation, Eulerian algorithms cover the volume of study with a fixed grid and compute the fluid transfer through the faces of the (fixed) grid cell volumes to follow the evolution of the system. Lagrangian formulations, on the other hand, compute the system by following the ever changing volume and shape of a particular individual element of gas (interestingly, the ‘Lagrangian’ formulation is also due to Euler eul1862 (1862), who employed this formalism in a letter to Lagrange, who later proposed these ideas in a publication by himself, lag1762 (1762); see Whitehurst Wh95 (1995)).

For a substantial part the success of the DTFE may be ascribed to the use of Delaunay tessellations as an optimally covering grid. This suggests that they may also be ideal for the use in moving grid implementations for hydrodynamical calculations. As in our SPH application, such hydrocodes with Delaunay tessellations at their core would warrant a close connection to the underlying matter distribution. Indeed, attempts towards such implementations have already been introduced in the context of a few specific, mainly two-dimensional, applications (Whitehurst Wh95 (1995), Braun & Sambridge BraSam95 (1995), Sukumar Suk98 (1998)). Alternative attempts towards the development of moving grid codes, in an astrophysical context, have shown their potential (Gnedin Gne95 (1995), Pen Pen98 (1998)).

For a variety of astrophysical problems it is indeed essential to have such advanced codes at one’s disposal. An example of high current interest may offer a good illustration. Such an example is the reionization of the intergalactic medium by the ionizing radiation emitted by the first generation of stars, (proto)galaxies and/or active galactic nuclei. These radiation sources will form in the densest regions of the universe. To be able to resolve these in sufficient detail, it is crucial that the code is able to focus in onto these densest spots. Their emphasis on mass resolution makes Lagrangian codes – including SPH – usually better equipped to do so, be it not yet optimally. On the other hand, it is in the low density regions that most radiation is absorbed at first. In the early stages the reionization process is therefore restricted to the huge underdense fraction of space. Simulation codes should therefore properly represent and resolve the gas density distribution within these voidlike regions. The uniform spatial resolution of the Eulerian codes is better suited to accomplish this. Ideally, however, a simulation code should be able to combine the virtues of both approaches, yielding optimal mass resolution in the high density source regions and a proper coverage of the large underdense regions. Moving grid methods, of which Delaunay tessellation based ones will be a natural example, may indeed be the best alternative, as the reionization simulations by Gnedin (Gne95 (1995)) appear to indicate. There have been many efforts in the context of Eulerian codes towards the development of Adaptive Mesh Refinement( AMR) algorithms( Berger BC89 (1989)), which have achieved a degree of maturity. Their chief advantage is their ability to concentrate computational effort on regions based on arbitrary refinement criteria, where, in the basic form at least, moving grid methods refine on a mass resolution criterion. However they are still constrained by the use of regular grids, which may introduce artifacts due to the presence of preferred directions in the grid. The advantages of a moving grid fluid dynamics code based on Delaunay tessellations have been most explicitly demonstrated by the implementation of a two-dimensional lagrangian hydrocode (FLAME) by Whitehurst (Wh95 (1995)). These advantages will in principle apply to any such algorithm, in particular also for three-dimensional implementations (of which we are currently unaware). Whitehurst (Wh95 (1995)) enumerated various potential benefits in comparison with conventional SPH codes, most importantly the following:

1. SPH needs a smoothing length h.
2. SPH needs an arbitrary kernel function W.
3. The moving grid method does not need an (unphysical) artificial viscosity to stabilize solutions.

The validity of the first two claims has of course also been demonstrated in this study for particle methods based on DTFE. Whitehurst showed additionally that there is an advantage of moving grid methods over Eulerian grid-based ones. The implementation of Whitehurst, which used a first-order solver and a limit on the shape of grid cells to control the effects of shearing of the grid, was far superior to all tested first-order Eulerian codes, and superior to many second-order ones as well. The adaptive nature of the Langrangian method and the fact that the resulting grid has no preferred directions are key factors in determining the performance of moving grid methods such as FLAME. For additional convincing arguments, including the other claims, we may refer the reader to the truly impressive case studies presented by Whitehurst (Wh95 (1995)).

## 6 Summary & discussion

Here we have introduced the DTFE as an alternative density estimator for particle fluid dynamics. Its principle asset is that it is fully self-adaptive, resulting in a density field reconstruction which closely reproduces, usually in meticulous detail, the characteristics of the spatial particle distribution. It may do so because of its complete independence of arbitrary user-specified smoothing functions and parameters. Unlike conventional methods, such as the kernel estimators used in SPH, it manages to faithfully reproduce the anisotropies in the local particle distribution. It therefore automatically reflects the genuine geometry and shape of the structures present in the underlying density field. This is in marked contrast with kernel based methods, which almost without exception produce distorted shapes of density features, the result of the convolution of the real structure with the intrinsic shape of the smoothing function. Its adaptive and local nature also makes it optimally suited for reconstructing the hierarchy of scales present in the density distribution. In kernel based methods the internal structural richness of density features is usually suppressed on scales below that of the characteristic (local) kernel scale. DTFE, however, is solely based upon the particle distribution itself and follows the density field wherever the discrete representation by the particle distribution allows it to do so. Its capacity to resolve structures over a large dynamic range may prove to be highly beneficial in many astrophysical circumstances, quite often involving environments in which we encounter a hierarchical embedding of small-scale structures within more extended ones.

In this study we have investigated the performance of the DTFE density estimator in the context of a Smooth Particle Hydrodynamics simulation of a multiphase interstellar medium of neutral gas. The limited spatial resolution of current particle hydrodynamics codes are known to implicate considerable problems near regions with e.g. steep density and temperature gradients. In particular their handling of shocks forms a source of considerable concern. SPH often fails in and around these regions, so often playing a critical and vital role in the evolution of a physical system. Our study consists of a comparison and confrontation of the conventional SPH kernel based density estimation procedure with the corresponding DTFE density field reconstruction method.

The comparison of the density field reconstructions demonstrated convincingly the considerable improvement embodied by the DTFE procedure. This is in particular true at locations and under conditions where SPH appears to fail. Filamentary and sheetlike structures provide telling examples of the superior DTFE handling with respect to the regular SPH method, with the most pronounced improvement occurring in the direction of the steepest density gradient.

Having shown the success of the DTFE, we are convinced that its application towards the analysis of the outcome of SPH simulations will prove to be highly beneficial. This may be underlined by considering a fitting illustration. Simulations of the settling and evolution of the X-ray emitting hot intracluster gas in forming clusters of galaxies do represent an important and cosmologically relevant example (see Borgani & Guzzo BorGuz01 (2001) and Rosati et al. Rosetal02 (2002) for recent reviews). The X-ray luminosity is strongly dependent upon the density of the gas. The poor accuracy of the density determination in regular SPH calculations therefore yields deficient X-ray luminosity estimates (see Bertschinger Ber98 (1998) and Rosati et al. Rosetal02 (2002) for relevant recent reviews). Despite a number of suggested remedies, such as separating particles according to their temperature, their ad-hoc nature does not evoke a strong sense of confidence in the results. Numerical limitations will of course always imply a degree of artificial smoothing, but by invoking tools based upon the DTFE technique there is at least a guarantee of an optimal retrieval of information contained in the data.

Despite its promise for the use in a variety of analysis tools for discrete data samples, such as particle distributions in computer simulations or galaxy catalogues in an observational context, its potential would be most optimally exploited by building it into genuine new fluid dynamics codes. Some specific (two-dimensional) examples of succesful attempts in other scientific fields were mentioned, and we argue for a similar strategy in astrophysics. One path may be the upgrade of current particle hydrodynamics codes by inserting DTFE technology. In this study, we have outlined the development of such a SPH-like hydrodynamics scheme in which the regular kernel estimates are replaced by DTFE estimates. One could interpret this in terms of the replacement of the user-specified kernel by the self-adaptive contiguous Delaunay cell, solely dependent on the local particle configuration. An additional benefit will be that on the basis of the localized connections in a Delaunay tessellations it will be possible to define a more physically motivated artificial viscosity term.

The ultimate hydrodynamics algorithm would combine the virtues of Eulerian and Lagrangian techniques. Considering the positive experiences with DTFE, it appears to be worthwhile within the context of ‘moving grid’ or ‘Lagrangian grid’ methods to investigate the use of Delaunay tessellations for solving the Euler equations. With respect to a particle hydrodynamics code, the self-adaptive virtues of DTFE and its ability to handle arbitrary density jumps with only one intermediate point may lead to significant improvements in the resolution and shock handling properties. Yet, for grid based methods major complications may be expected in dealing with the non-regular nature of the corresponding cells, complicating the handling of flux transport along the boundaries of the Delaunay tetrahedra.

The computational cost of DTFE resembling techniques is not overriding. The CPU time necessary for generating the Delaunay tessellation corresponding to a point set of $`N`$ particles is in the order of $`O(N\mathrm{log}N)`$, comparable to the cost of generating the neighbour list in SPH. Within an evolving point distribution these tessellation construction procedures may be made far more efficient, as small steps in the development in the system will induce a correspondingly small number of tetrahedron (identity) changes. Such dynamic upgrading routines are presently under development.

In summary, in this work we have argued for and demonstrated the potential and promise of a natural computational technique which is based upon one of the most fundamental and natural tilings of space, the Delaunay tessellation. Although the practical implementation will undoubtedly encounter a variety of complications, dependent upon the physical setting and scope of the code, the final benefit of a natural moving grid hydrodynamics code for a large number of astrophysical issues may not only represent a large progress in a computational sense. Its major significance may be found in its ability to address fundamental astrophysical problems in a new and truely natural way, leading to important new insights in the workings of the cosmos.

###### Acknowledgements.

The authors thank Vincent Icke for providing the incentive for this project and for fruitful discussions and suggestions. WS is grateful to Emilio Romano-Díaz for stimulating discussions. RvdW thanks Dick Bond and Bernard Jones for consistent inspiration and encouragement to investigate astrohydro along tessellated paths. We are also indebted to Jeroen Gerritsen and Roelof Bottema for providing software and analysis tools.
This work was sponsored by the stichting Nationale Computerfaciliteiten (National Computing Facilities Foundation) for the use of supercomputer facilities, with financial support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Netherlands Organisation for Scientific Research, NWO)


format: text

text:
# Improved Tight Binding Parametrization for the Simulation of Stacking Faults in Aluminium

## Abstract

We refit the NRL tight binding parameterization for Aluminium by Mehl *et al* \[Phys. Rev. B, 61, 4894 (2000)\], to a database generated via full potential Linearized Augmented Plane Wave (LAPW) Density Functional Theory (DFT) calculations. This is performed using a global optimization algorithm paying particular attention to reproducing the correct order of the angular symmetries of the tight binding fcc and bcc bandstructure. The resulting parameterization is found to better predict the hcp phase and both the stable and unstable planar stacking fault defect energies.

The empirical tight binding technique can be seen as a compromise between the very accurate but computationally expensive *ab initio* Density Functional Theory (DFT) methods, and the fast but less accurate empirical potential methods such as the Embedded Atom Method (EAM). Although great reductions in computational cost has been achieved for *ab initio* methods by employing a minimal basis set and efficient algorithms that scale linearly with atom number Solér et al. (2001), fully self-consistent calculations are still relatively slow for “large” scale Molecular Dynamics (MD) simulations. For empirical tight binding methods, a range of linear scaling techniques have also been developed Bowler et al. (1997) and the development of accurate tight binding parameterizations is of utmost importance, if this method is to remain competitive. Furthermore, the strength of the tight binding method, contrary to many EAM parameterizations, is its ability to accurately predict a range of structural properties not included in the fitting database.

The NRL tight binding formalism is a total energy tight binding method based on the fitting of bandstructure eigenvalues and total energies from *ab initio* calculations. One of the advantages of this method is that one is able to fit to the bandstructure and total energy of a given structure simultaneously be redefining the energy zero of the sum of eigenvalues. From DFT, the total energy can be expressed as the sum of the single electron eigenvalues plus an ionic-like term

$$E(n(𝐫))=\underset{i}{}f(\mu ϵ_i)ϵ_i+F(n(𝐫))$$
(1)
where $`f(x)`$ is the Fermi function, $`ϵ_i`$ the electronic eigenvalues, $`\mu `$ the chemical potential, and $`F(n(𝐫))`$ the so-called “correction term” Mehl and Papaconstantopoulos (1996); Sang H Yang and Papaconstantopoulos (1998). Following Mehl *et al* this can be reformulated by shifting the energy zero of the bandstructure term by a value exactly canceling the “correction term”:

$$E(n(𝐫))=\underset{i}{}f(\mu ^{}ϵ_i^{})ϵ_i^{}$$
(2)
This shift is taken into account by making the onsite terms environmentally dependent. Using a nonorthogonal two-center Slater-Koster formulation Slater and Koster (1954), the tight binding model is then parameterized by fitting the bandstructure and total energies as a function of lattice constant and crystal structure of the tight binding model to a database produced by *ab inito* calculations. The philosophy of the NRL formalism is to include only a handful of high symmetry structures in the fitting database. This method has proven to be successful for a wide range of elements. Further details can be found in Refs. Mehl and Papaconstantopoulos (1996, 1998).

We have refitted a tight binding parameterization by Mehl *et al* Mehl et al. (2000a) for Aluminium using an $`s`$, $`p`$ and $`d`$ orbital basis, to the full potential LAPW DFT calculations of WIEN2kBlaha et al. (1999). As an alternative to the conventional least squares procedure we have used Adaptive Simulated Annealing (ASA), a stochastic global optimization algorithm Frøseth et al. (2001); Ingber (1993). Assuming a probability density of states $`g_T(𝐱)`$ for the parameters space $`\{x_i\}`$, and a cost function defined by the “energy” $`E(x)`$. At any step, $`k`$, the “probability” for accepting the new cost function $`E_{k+1}`$ is given by:

$$h(\mathrm{\Delta }E)=\{\begin{array}{cc}\mathrm{exp}(\mathrm{\Delta }E/T)\hfill & \text{for }\mathrm{\Delta }E>0,\hfill \\ 1\hfill & \text{for }\mathrm{\Delta }E<0.\hfill \end{array}$$
(3)
where the parameter $`T=T(k)`$ defines the “annealing schedule”. For such a system the principle of detailed balance holds, which means that, in theory, all states of the system will be sampled by the algorithm and a global minimum can be found. A weakness of the simulated annealing approach is that it needs to sample a large part of parameter space and can therefore be inherently slow. This is partly remedied in the ASA algorithm by assuming a probability density for the parameters with “fat” logarithmic tails

$$g_T(x)\underset{i}{}g_T^i(x_i)=\underset{i}{}\frac{1}{2(|x_i|+T_i)\mathrm{ln}(1+T_i)}$$
(4)
where $`x=\{x_i[1,1];i=1,N\}`$ is the normalized parameter set. This allows for a very fast annealing schedule

$$T_i(k)=T_{t0}\mathrm{exp}c_ik^{1/N}$$
(5)
in contrast to the conventional linearly decreasing annealing schedule. $`T_{i0}`$ and $`c_i`$ are free parameters, which allows the algorithm to adapt to a changing “environment” for the cost functions as the search proceeds. Further details can be found in ref. Frøseth et al., 2001.

In the fitting database we included band-structure data as a function of lattice constant for the the fcc and bcc crystal structures using irreducible k-point meshes of 146 for fcc and 91 points for bcc. In order to ensure the correct angular character for each eigenvalue and further constrain the Slater-Koster parameters we block-diagonalized the Hamiltonian for a large number of high symmetry points and directions: the $`\mathrm{\Gamma }`$, $`X`$, $`L`$, $`W`$ points and $`\mathrm{\Lambda }`$, $`\mathrm{\Delta }`$, $`\mathrm{\Sigma }`$, $`S`$, $`Z`$ directions for the fcc structure, and the $`\mathrm{\Gamma }`$, $`H`$, $`P`$, $`N`$ points and $`\mathrm{\Lambda }`$, $`\mathrm{\Delta }`$, $`\mathrm{\Sigma }`$, $`F`$, $`D`$, $`G`$ directions for the bcc structure.

For the DFT calculations we used the input parameters: $`R_{mt}=2.5`$ Bohr, $`R_{mt}K_{max}=7`$ and $`G_{max}=14`$ $`\mathrm{Ry}^{1/2}`$. Here, $`R_{mt}`$ is the muffin tin radius, $`K_{max}`$ is the plane wave cut-off and $`G_{max}`$ is the maximum Fourier component of the electron density. For the exchange-correlation potential we used the Generalized Gradient Approximation (GGA) of Perdew et. al Perdew et al. (1996).

Fig. 1 displays the TB bandstructure for fcc at $`a=7.65`$ Bohr compared to FLAPW *ab inito* calculations and the TB parameterization of Mehl *et. al*. In contrast to the previous parametrization, the ordering with respect to angular character of all bands below 1.25 Ry is now correct. This is important, since otherwise the relative contribution to the bandstructure from each atomic orbital in the LCAO-expansion would break the basic symmetry-requirements of the crystal potential Inui et al. (1995).

Fig. 2 displays the energies of state for the fcc, hcp and bcc structures as obtained from a second order Birch fit to the data points derived from the present TB model. Also shown are FLAPW *ab inito* data for fcc and bcc. Tab. 1 compares the elastic constants from the TB parameterization of the present work with the parameterization of Mehl *et al*Sang H Yang and Papaconstantopoulos (1998), *ab initio* DFT calculations using WIEN2k and experimental values. As can be seen, the TB values are in excellent agreement with the *ab inito* and experimental values, where only the bulk modulus is overestimated. This later discrepancy is probably due to more weight being given (in the present fitting process) to bandstructure data than the equation of state total energy data. The hcp data is in effect a prediction of the tight binding parameterization since only the fcc and bcc structures where fitted. At the obtained equilibrium lattice constants, our fit predicts a fcc-hcp energy difference of 2.1 mRyd per atom, which agrees well with that of the converged (with respect to k-point mesh density using the previously stated *ab initio* parameters) WIEN2k value of 2.0 mRyd. The TB parametrization of Mehl *et. al.* returns a value of 1.8 mRyd.

An important fcc material property for any interatomic model to reproduce is the intrinsic stacking fault surface energy density since it plays a crucial role in the disassociation width of a full dislocation into its constituent Shockley partials and is believed to play a role in the issue of why in some cases only partial dislocations are observed in atomistic modeling of the mechanical properties of nanocrystalline materials Yamakov et al. (2001); Swygenhoven et al. (2002). An unrelaxed stacking fault is modeled by stacking layers of closed packed planes in the  direction using a particular stacking sequence. We choose the primitive vectors of the unit cell to be:

$$\begin{array}{cc}\hfill 𝐚_\mathrm{𝟏}=& \frac{1}{2}a𝐲+\frac{1}{2}a𝐳\hfill \\ \hfill 𝐚_\mathrm{𝟐}=& \frac{1}{2}a𝐱+\frac{1}{2}a𝐳\hfill \\ \hfill 𝐚_\mathrm{𝟑}=& (4+\frac{q}{6})a𝐱+(4+\frac{q}{6})a𝐲(4\frac{q}{3})a𝐳\hfill \end{array}$$
(6)
where $`a`$ is the lattice constant. $`q`$ is a variable controlling the slip in the  direction. This is equivalent to a supercell of 12 close-packed planes, with a variable slip at the 12th layer due to the periodic boundary conditions. For a $`q=0`$ we get a perfect fcc crystal represented by the stacking sequence $`ABCABC`$, and for $`q=1`$ we get the full stacking fault $`ABC|BC`$, where the stacking sequence at the interface corresponds to two 111 planes in which the atoms are hcp coordinated. A mapping of the structural energy for all values of $`q`$ from $`0`$ to $`1`$ is called a *Bain Path* (or also a generalized stacking fault curve) where an intrinsic staking fault occurs at $`q=1`$. Fig. 3 shows the stacking fault energy as a function of the slip variable, $`q`$ relative to the fcc energy at $`q=0`$. The maximum of the curve defines the unstable stacking fault configuration, which in fig. 3 is near $`q=0.6`$.

Table 2 presents the calculated values for the intrinsic stacking fault energy, $`\gamma _{is}`$, unstable stacking fault energy, $`\gamma _{us}`$, and the ratio $`\gamma _{is}/\gamma _{us}`$, compared to the TB parameterization of Mehl et. alMehl et al. (2000b). In addition we include published results from DFT and EAM calculations, and experimental values. The experimentally derived values for $`\gamma _{is}`$ vary from a low value of 110 $`\mathrm{mJ}/\mathrm{m}^2`$ to a high value of 280 $`\mathrm{mJ}/\mathrm{m}^2`$. However, *ab initio* DFT calculations have consistently shown results in the range 161-166 $`\mathrm{mJ}/\mathrm{m}^2`$. The EAM potential of Voter and Chen shows values for the stacking fault energies that are seriously underestimated. The important point here is that these values were not included in the fitting database for this potential Voter and Chen (1987). Similar results have been produced by the EAM potential of Ercolessi and Adams, which has been used in large scale Molecular Dynamics (MD) simulations of extended dislocations Yamakov et al. (2001). In contrast to this, the recently published EAM potential of Mishin and Farkas, show values for the intrinsic and unstable stacking fault energies which are more in line with the *ab initio* calculations. The reason for this is that a value for the intrinsic stacking fault energy was included in the fitting database Mishin et al. (1999). The TB parameterization of the present work shows a marked improvement over the older parameterization of Mehl *et al* Especially the intrinsic stacking fault energy, $`\gamma _{is}`$, with a value of 146 $`\mathrm{mJ}/\mathrm{m}^2`$ which is in much better agreement with *ab inito* calculations than a value of 97 $`\mathrm{mJ}/\mathrm{m}^2`$ of the original fit. Another important improvement is the value of the ratio $`\gamma _{us}/\gamma _{is}`$, which is a measure of the relative cost of dislocation nucleation with respect to the intrinsic stacking fault energy. For the new TB parameterization this value is 1.3 which is in good agreement with the value of 1.4 from *ab inito* calculations.

The improved prediction of the stacking fault energies is consistent with the corresponding improvement of the hcp cohesive energy and lattice constant of the present model. At the nearest neighbor bond length, the local fcc and hcp environments differ only with respect to their dihedral angle arrangements. In terms of energy and lattice constant this difference should be sensitive to the angular character of, particularly, the $`\pi `$ and $`\delta `$ bonds of the $`p`$ and $`d`$ states, and it is precisely such bonding elements that the current parameter-set more accurately reproduces at and near the Fermi level.

In conclusion, we have refitted the tight binding parameters of Mehl *et al* for Aluminium to *ab initio* FLAPW eigenvalues. We have used a stochastic global optimization algorithm and an extended database of high symmetry eigenvalues. As a result we obtain a parameterization with an improved representation of the angular symmetries for the Slater Koster parameters. As a consequence we obtain improved predictions for the stacking fault energies without including this information explicitly in the fitting database.

Parts of this work has been supported by the Norwegian Research Council through CPU time on the NOTUR Origin 3800. The static tight-binding code was provided by M. J. Mehl of the Naval Research Laboratory under the U.S. Department of Defense CHSSI program.


text:
# Surface-active dust in a plasma sheath

## Abstract

The inhomogeneity of the plasma pressure near a conducting electrode is a cause for introducing the surface tension associated with the plasma-electrode interface. We evaluate the dependence of the surface tension on the density of the charged dust immersed in the plasma sheath. In a wide range of parameters, the surface tension turns out to be an increasing function of the dust density.

Throughout the recent decade, an admixture of charged macro-particles, electrons, ions, and neutral atoms, called a dusty plasma, has been the subject to some thousands of studies. Numerous aspects of the physics of dusty plasmas are discussed in the series of reviews starting with TMT02  and the monograph shukla02 .

Many laboratory experiments deal with relatively small number of dust grains levitating above a horizontal rf-powered or dc-biased electrode. Under these conditions, dust grains are negatively charged and the gravity force is compensated by the strong electric field in a plasma sheath adjacent to an electrode. Dust suspended in a plasma sheath self-organizes itself in various single or multi-layered structures. In particular, recent experiments tue  demonstrated that with growing rf power or gas pressure there evolves a void, that is, a dust-free region, in the center of a single dust layer. Since negatively charged grains levitating at the same height repulse one another via the screened Coulomb potential, there are no evident reasons for the void emergence. The explanation of this phenomenon given in tue  implies existence of an attractive interaction between grains; however, the nature of the interaction remains obscure.

To avoid confusion, it should be stressed that three-dimensional voids, which have been observed earlier goree ; therm ; mg1 ; PKE , evolve in the plasma bulk. The theory of three-dimensional voids theory1 ; theory2  demands large number of grains and strong influence of the dust upon the discharge structure. At first sight, both assumptions fail under the conditions of the experiment tue , where some hundreds of grains only have been used.

The main purpose of this Letter is to introduce the concept of plasma surface tension. We demonstrate that even a small number of dust grains can modify the plasma sheath profile in such a way that the sheath may confine the dust in the horizontal direction. Our reasoning is based on the analysis of the plasma pressure in the sheath-presheath area of a gas discharge. We show that the horizontal electric force appearing due to the inhomogeneity of a dust layer is conveniently treated in terms of the plasma surface tension. It should be stressed that explanation of the experiment tue  is beside our purpose: there is evidently a number of relevant factors that are not discussed here.

The schematic of the model implemented here is depicted in Fig. 1. The plasma consisting of electrons and singly charged ions is described with the help the continuity equation

$$\frac{dn_i(z)v_i(z)}{dz}=\nu _{ion}n_e(z),$$
(1)
the momentum balance equation

$$\frac{dm_in_i(z)v_i^2(z)}{dz}+en_i(z)\frac{d\varphi (z)}{dz}=0,$$
(2)
and Poisson’s equation

$$\frac{d^2\varphi (z)}{dz^2}=4\pi e(n_e(z)n_i(z))4\pi \rho _d(z).$$
(3)
Here, $`\nu _{ion}`$ is the electron impact ionization frequency and $`\rho _d(z)`$ is the charge density of aerosol grains. Other notation is standard. Electrons are supposed to be Boltzmannian: $`n_e(z)=n_0\mathrm{exp}\left(e\varphi (z)/T_e\right)`$, where $`n_0`$ is a normalizing constant. We also assume that the mean free path of the ion-neutral collisions is large compared to the Debye length, which determines the spatial scale of the sheath.

Since the electron and ion fluxes at the electrode ($`z=0`$) are equal, one of the boundary conditions to Eqs. (1-3) is

$$n_i(0)v_i(0)+n_e(0)\sqrt{\frac{T_e}{2\pi m_e}}=0.$$
(4)
Far from the electrode the plasma is quasi-neutral, i.e., $`n_e=n_i=n_0`$. Generally, the equilibrium plasma density, $`n_0`$, is the eigenvalue of the problem and it is determined by the balance of ionization and losses of charged particles. For this reason, the bulk plasma density may depend on the dust charge located in the sheath because the latter modifies the electric field profile, which, in its turn, influences the net particle flux at the electrode. However, here we assume the overall size of the plasma to be large enough, so that the equilibrium density, $`n_0`$, is conditioned by the recombination in the plasma volume. Therefore, the plasma density tends to a certain fixed value at $`z\mathrm{}`$. In order to patch the quasi-neutral plasma and the sheath, the equations (1,2) are solved assuming that $`n_e(z)=n_i(z)`$ and $`\rho _d(z)=0`$ (e.g., riemann ) that results in

$$n_i(z)=\frac{n_0v_s^2}{v_s^2+v_i(z)^2},$$
(5)
where $`v_s=\sqrt{T_e/m_i}`$ is the ion sound velocity. It is assumed that at the end of the simulation area, $`z=z_0`$, the ion flow velocity is fixed, i.e., $`v_i(z_0)=v_0`$. The value of $`v_0`$ should be small enough; all solutions exemplified below are obtained for $`v_0=0.01v_s`$. Then, Eqs. (1,5) are used to evaluate boundary values of the density, $`n_i(z_0)`$, the potential, $`\varphi (z_0)`$, and its derivative, $`\varphi ^{}(z_0)`$; this results in four boundary conditions for Eqs. (1-3). The length of the simulation area, $`z_0`$, is determined with the help of the fifth boundary condition (4).

In the framework of our model, the dust is represented as an infinitesimally thin charged massive layer levitating at $`z=z_l`$, that is, its charge density is

$$\rho _d(z)=q\sigma \delta (zz_l),$$
(6)
where $`q`$ ($`q>0`$) is the charge of a single grain and $`\sigma `$ is the surface dust density. The equilibrium levitation height, $`z_l`$, depends on the vertical force balance. The largest forces acting upon dust in a plasma sheath are the electric field force and the gravity force. Since the $`z`$-component of the electric field, $`E(z)=\varphi ^{}(z)`$, is discontinuous at $`z=z_l`$, the vertical force balance is written as

$$\frac{1}{2}q\left[E(z_l+0)+E(z_l0)\right]=Mg,$$
(7)
where $`M`$ is the mass of a single grain and $`g`$ is the gravity acceleration (in Fig. 1, the Earth is at the bottom). On the other hand, the electric field discontinuity is given by

$$E(z_l+0)E(z_l0)=4\pi q\sigma $$
(8)
and the potential, $`\varphi (z)`$, is a continuous function at $`z=z_l`$. The relations (7,8) allows us to patch the numeric solutions of Eqs. (1-3) below ($`0<z<z_l`$) and above ($`z_l<z<z_0`$) the layer, and to obtain the equilibrium levitation height, $`z_l`$, and the net thickness of the sheath-presheath area, $`z_0`$.

It is well-known that the grain charge, $`q`$, in Eqs. (7,8) depends on the plasma parameters, and its variability should be also taken into account vc00 . However, here we assume that the grain charge is fixed; the results of preliminary computer runs with the variable grain charge did not alter our main conclusions. We also ignore the ion wind force in Eq. (7) as long as absorption of plasma by dust grains. Simple estimations show that these processes are of minor importance.

Our main objective is to analyze the plasma pressure profile in the sheath-presheath area. Evidently, Eqs. (1-3) ensure the conservation of the net linear momentum of the system. The tensor of the momentum flux, or the pressure tensor is a sum of the ion momentum flux, the electron pressure and Maxwell’s stress: $`P_{\alpha \beta }=m_in_iv_{i\alpha }v_{i\beta }+n_eT_e\delta _{\alpha \beta }+𝑬^2\delta _{\alpha \beta }/8\pi E_\alpha E_\beta /4\pi `$ ($`\alpha ,\beta =x,y,z`$). In the one-dimensional model described above, the pressure tensor is diagonal: $`P_{zz}=P_l`$ and $`P_{xx}=P_{yy}=P_{tr}`$. The longitudinal pressure is

$$P_l(z)=m_in_i(z)v_i(z)^2+n_e(z)T_e\frac{1}{8\pi }E(z)^2.$$
(9)
Since within our model there are no momentum losses in the plasma, $`P_l`$ is constant in the absence of dust and it is equal to the plasma electron pressure at $`z\mathrm{}`$, $`P_l=n_0T_e`$.

The transverse pressure is written as

$$P_{tr}(z)=n_e(z)T_e+\frac{1}{8\pi }E(z)^2.$$
(10)
In contrast with Eq. (9), there is no contribution of the ion momentum flux and the electric field pressure is of the opposite sign.

Although the transverse pressure tends to the electron pressure in the plasma bulk, $`P_{tr}|_z\mathrm{}=n_0T_e`$, it is essentially nonuniform in the sheath-presheath area adjacent to the electrode. This circumstance allows us to define the surface tension at the plasma-electrode interface as

$$s(\sigma )=\underset{0}{\overset{\mathrm{}}{}}𝑑z\left(n_0T_eP_{tr}(z)\right).$$
(11)
The physical meaning of this quantity is fairly evident. Since the bulk plasma pressure is positive, two parts of the discharge repulse one another. As we shall demonstrate later, the transverse pressure (10) is reduced near the wall; therefore, there arises the additional horizontal force attracting two parts of the plasma-wall interface. The surface tension (11) characterizes the magnitude of this force per unit length. It should be noted that usually surface tension is defined in terms of surface free energy. The latter is obviously meaningless in application to an open, nonequilibrium system like the plasma-wall boundary. However, due to the momentum conservation, the surface tension given by (11) is a well-defined concept. Of particular interest for our purposes is the dependence of the surface tension (11) on the dust density.

The utility of the notion (11) may be illustrated in the following way. Suppose that the layer density is a weakly varying function of the transverse coordinate, $`x`$, that is, $`\sigma /xϵ`$, where $`ϵ1`$. Integrating $`P_{xx}`$ component of the pressure tensor with respect to $`z`$ at nearby points $`x`$ and $`x+\mathrm{\Delta }x`$ we see that the horizontal force (i.e., the momentum flux) acting upon a small part of the layer of the width $`\mathrm{\Delta }x`$ is written in two identical forms:

$$F_x=q\sigma E_x\mathrm{\Delta }x\mathrm{\Delta }x\frac{}{x}𝑑zP_{xx}(x,z).$$
(12)
Comparing both parts of the identity and noting that $`E_x/xO(ϵ)`$ we conclude that in order to obtain the first-order term of the expansion of the horizontal force in powers of $`ϵ`$ it is sufficient to evaluate the surface tension (11) ignoring the spatial inhomogeneity of the layer density. Therefore, $`f_x=s(\sigma (x))/x`$. This simple reasoning may be confirmed by the more tedious asymptotic expansion of the hydrodynamic equations in powers of $`ϵ`$.

We can also extract some useful information about the horizontal force even for the sharply bounded layer. Suppose that the dust layer occupies the area $`x<0`$, as depicted in Fig. 1. The levitation height, $`z_l`$, now depends on $`x`$, the electric field and the ion flow near the edge of the layer are essentially two-dimensional. However, the horizontal component of the force is the difference $`f_x=s_0s_1`$, where $`s_{0,1}`$ are given by the integral (11) evaluated across the sheath in the dust-free area, at $`x=x_0`$, and in the presence of dust, $`x=x_1`$ (Fig. 1). Assuming that both points $`x_{0,1}`$ are far enough from the layer edge, the horizontal components of the ion velocity and the electric field are negligible. Therefore, the values of the surface tension, $`s_{0,1}`$, which determine the net force exerted upon the edge, are provided by the solution of the one-dimensional problem (1-3).

Before discussing the results of the numeric evaluation of the surface tension, it is convenient to introduce the normalized quantities. In the following, the coordinate, $`z`$, is measured in the units of the electron Debye length, $`\lambda _D=\sqrt{T_e/4\pi e^2n_0}`$, the pressure is normalized to $`n_0T_e`$, and the surface tension is normalized to $`\lambda _Dn_0T_e`$. The dimensionless surface density of the dust layer is $`\eta =q\sigma /(2en_0\lambda _D)`$ and the normalized grain mass is $`\mu =Mge\lambda _D/(qT_e)`$. Under the typical experimental conditions, $`\mu 1`$, while $`\eta 0.1`$. All computations discussed below have been performed with the sufficiently small value of the ionization frequency, $`\nu _{ion}/\omega _{pi}=0.01`$, where $`\omega _{pi}`$ is the ion plasma frequency.

Examples of the pressure and the electric field profiles are depicted in Fig. 2. The net simulation length for this particular run is $`z_0=68.5`$. In order to make the effect discernible in the figure, we use here a sufficiently large value of the normalized dust density, $`\eta =0.1`$. The longitudinal pressure drop (the curve 1 in Fig. 2a) at the charged layer is compensated by its weight; outside of the layer, the longitudinal pressure is constant. Contrary, the transverse pressure (10) is strongly inhomogeneous inside the sheath (the curves 2 and 3 in Fig. 2a). Of importance is that the dust reduces the transverse pressure between the layer and the electrode. As it is readily seen from Fig. 2b, the effect is mostly due to the reduction of the electric field.

The vertical dashed lines in Fig. 2 show the position of the sheath-presheath interface, $`z_s`$, for $`\eta 0`$, which is conditionally defined as the point where the ion flow velocity is equal to the ion sound velocity. The dependence of the coordinate $`z_s`$ on the dust parameters is rather weak. For this particular example, the dust layer levitates in the supersonic ion flow, $`z_l<z_s`$; with reducing grain weight, $`\mu `$, it moves up to the subsonic area, $`z_l>z_s`$.

Fig. 3 shows the dependence of the surface tension on the dust layer density, $`\eta `$, for various values of the grain mass, $`\mu `$. The curves 1-3 in Fig. 3 correspond to the dust layer levitating in the supersonic flow, $`z_l<z_s`$, the curve 4 corresponds to lighter dust grains situated in the presheath area, $`z_l>z_s`$. The depicted curves show that if the dust layer is located in the sheath, the surface tension increases with $`\eta `$, while with the dust layer shifted to the presheath, the surface tension may become a decreasing function of $`\eta `$. This is also shown in Fig. 4, where the dependence of the derivative $`s^{}(\eta )|_{\eta =0}`$ on the grain mass, $`\mu `$, is plotted.

As we have already mentioned, the dependence of the surface tension on the dust density determines the force acting upon the edge of the semi-infinite layer. If $`s(\eta )>s(0)`$ then the force $`f_x=s(0)s(\eta )`$ is negative, as shown in Fig. 1, that is, the surface tension at the plasma-wall interface tends to expel the dust. Otherwise, the dust is drawn in the sheath. As it is seen from Figs. 3,4, in dependence on the grain mass and the dust density both alternatives are possible.

Here, we treat dust layer as a rigid incompressible medium. However, the short-range interaction between the dust grains result in the surface pressure, $`p_s(\sigma )`$, which also depends on the density hebner . The surface tension (11) reduces the surface pressure, that is, the net equation of state of the dust layer is given by $`p=p_s(\sigma )s(\sigma )`$. Since $`p_s(\sigma )`$ exponentially tends to zero with decreasing $`\sigma `$, the compressibility, $`dp/d\sigma `$, of a rarefied dust layer is layer is negative if $`s^{}(0)>0`$. Such a medium is unstable. The detailed investigation of this instability and arising structures will be discussed elsewhere.

To summarize, charged dust may be treated as a surface-active substance that is capable of altering the surface tension at the plasma-electrode interface. Unlike usual surfactants, the dust increases the surface tension in a wide range of parameters. Although we did not discuss here the emergence of the two-dimensional void tue , it should be pointed out that a similar effect is well-known in liquids. To observe a void, it is enough to drop some liquid soap at the surface of water with floating fine powder.

This study was supported in part by the Netherlands Organization for Scientific Research (grant no. NWO 047.008.013) and the Russian Foundation for Basic Research (project no. 02-02-16439).