text:
# Vector- and tensor-valued descriptors for spatial patterns

## 0.1 Introduction

The spatial patterns originating from the polymorphic aggregation of matter in Nature occur on vastly different length scales and with unlimited variety. For instance, the pattern shown in Figure 1 could represent a biological structure obtained with X-ray tomography, or a colloidal cluster in a milky emulsion. In fact, Figure 1 displays an isodensity contour of cosmic matter as traced by the galaxies surrounding the Milky way, which lies in the center of the plot.
The visual impression from such pictures may already convey valuable insight into the dynamic evolution or – in the context of biology – into the functionality of such structures. However, for unbiased inferences

as well as for the comparison with model simulations, one needs objective and quantitative measures to characterize the geometry and topology of typical structural motifs.
In this article we deal with a family of functions which provide a unique description of the morphology displayed by spatial patterns. We call these functions the Minkowski valuations; they are extensions of the well-known Minkowski functionals often mentioned within this volume (see the contributions by C. Arns et al. and M. Lösche and P. Krüger in this volume as well as beisbart\_mecke:additivity  in the previous volume LNP 554). Therefore they share some of the simplicity and beauty which the Minkowski functionals own. In particular, they are characterized by simple axioms which set intuitively reasonable standards for a morphological description. In the same way as the Minkowski functionals generalize the notion of the volume, the Minkowski valuations extend the notion of the mass vector (or center of mass) and of the inertia tensor. Like their scalar relatives the majority of the Minkowski valuations are curvature measures integrated over the surface of a pattern.
We start with a short motivation in Section 0.2 by showing that there are good reasons to go beyond the scalar Minkowski functionals. After having outlined ways to generalize the Minkowski functionals (Section 0.3), we review some of the extensions’ mathematical properties in Section 0.4. To illustrate how useful the Minkowski valuations are, we discuss concrete applications covering submolecular up to cosmological scales (Sections 0.5 and 0.6).
Before, two cautionary remarks are appropriate: first, in order to distinguish conceptually between the scalar Minkowski functionals and their higher-rank extensions, we shall address the whole entirety of Minkowski functionals and their extensions via the term “Minkowski valuations” (MVs), whereas “Minkowski functionals” (MFs) is reserved to the scalar Minkowski functionals as before. Secondly, there is a vivid theoretical development of higher-rank Minkowski valuations these days beisbart\_schneider:tensor ; beisbart\_schneider:tensor2 ; beisbart\_beisbart:tensor . Therefore, this overview can not provide a complete survey on Minkowski valuations; on the contrary, substantial progress is to be expected both from the mathematical side as well as from applications. For instance, morphological analysis based on Minkowski valuations may be useful to quantify molecular orientations (see the contributions by F. Schmid and N. H. Phuong, U. S. Schwarz and G. Gompper, D. Jeulin, C. Arns et al., and H.-J. Vogel in this volume). Our vector-valued descriptors may also be used to refine the morphological analysis of Langmuir monolayer phases (presented by M. Lösche and P. Krüger in this volume) which has so far been based merely on scalar MFs. Another application may cover the geometric characterization of liquid foams described by F. Graner in this volume, where stress tensors need to be related to a complex spatial structure.

## 0.2 Going beyond scalar Minkowski functionals – a physical motivation

The following problem is typical for challenges occurring in physics: the formation of galaxy clusters depends in a complicated way on the background cosmology usually described in terms of a cosmological model. Such a model is defined through the values of the cosmological parameters and the power spectrum of the initial density perturbations beisbart\_peebles:principles ; beisbart\_bartelmann:arcIV . In order to investigate the impact of the background cosmology on galaxy clusters one carries out gravitational $`N`$-body simulations; typical results comprise cluster images such as those shown in Figure 2. To compare the models and to check which model fits best present-day cluster observations, one has to use quantitative and discriminative morphological measures.

The Minkowski functionals certainly place practical tools at the physicist’s disposal. Defined on the convex ring $`𝒦_d`$, i.e. for finite unions of compact, convex sets in $`d`$-dimensional Euclidean space, they are the only motion-invariant, convex-continuous, and additive descriptors. Despite this very general characterization there are only a finite number of Minkowski functionals, namely $`(d+1)`$ in $`d`$ dimensions (characterization theorem, beisbart\_hadwiger:vorlesung ; beisbart\_klain:hadwiger ). They have simple geometrical meanings derived from their representation as integrals. These integrals either extend over the pattern itself to yield the volume, or cover its surface which is weighted with combinations of the local curvatures in order to give – in two dimensions, e.g. – the perimeter and the integrated curvature (which equals the Euler characteristic). The Minkowski functionals also arise quite naturally as expansion coefficients in the Steiner formula beisbart\_schneider:brunn  specifying the volume of the parallel body of a convex body. Indeed, the Minkowski functionals have been extensively applied to physical structures at vastly different length scales, see beisbart\_mecke:additivity  for an overview and beisbart\_kerscher:statistical  for astrophysical applications. On the mathematical side, numerous useful results such as the principal kinematic formulae or the Crofton formulae (see beisbart\_hadwiger:vorlesung ) have been proven with integral geometric methods.
However, the morphometry with Minkowski functionals is by no means comprehensive. Combining two of their defining principles, additivity and homogeneity, one can easily see, that, e.g., moving a connected part of a pattern results in significant visible changes which can not be detected by the Minkowski functionals at all. This is illustrated in the first panel of Figure 3 representing schematic clusters. It is surely of interest, however, in which way the subclumps of a system are arranged, and how they are aligned with respect to each other (see the second panel of Figure 3). Therefore, we have to move beyond the scalar Minkowski functionals. It turns out that one can generalize the Minkowski functionals without leaving the mathematical framework of integral geometry.

## 0.3 The hierarchy of Minkowski valuations – extending the framework

There are several ways to extend the framework of the Minkowski functionals. For instance, one may start with Steiner’s formula.

#### The Steiner formula.

It is a natural choice to employ the volume $`V(K)=_KdV`$ of a body $`K`$ for its very coarse description. If the parallel body $`K_ϵ`$ is considered, which contains all points closer to $`K`$ than $`ϵ>0`$, and if $`K`$ is convex, then one finds from Steiner’s formula an expansion for $`V(K_ϵ)`$ in powers of $`ϵ`$,

$$V\left(K_ϵ\right)=\underset{\nu =0}{\overset{d}{}}\left(\genfrac{}{}{0pt}{}{d}{\nu }\right)ϵ^\nu W_\nu ,$$
(1)
where the Minkowski functionals $`W_\nu `$ arise as expansion coefficients. The latter functionals are useful morphometric concepts and can be extended to non-convex patterns. In this sense, the volume quite generally leads to a family of related measures.
Since we saw that in some applications (recall the upper panels of Figure 3) the position of a body may be worth considering, we specify the location of $`K`$ with its mass vector

$$V^1\left(K\right)_KdV𝐱.$$
(2)
Alternatively, we could introduce its center of mass $`𝐩_0V^1/V`$.
Moreover, as the lower panels of Figure 3 indicate, the body’s orientation may be relevant. In this case we have to move on to even higher moments; for instance, the tensor

$$V^2\left(K\right)_KdV\mathrm{𝐱𝐱}$$
(3)
(where $`\mathrm{𝐱𝐱}`$ denotes the symmetric tensor product111As beisbart\_beisbart:tensor  show, there are no non-trivial antisymmetric Minkowski valuations of rank $`r1`$ at all, i.e. no motion-covariant, additive, and conditionally continuous tensor functions with rank $`r1`$, see below Section 0.3 for the explanation of these concepts.) contains information on main directions of the body $`K`$. Note, that both the mass vector of a body and its inertia tensor (the traceless part of Equation 3: $`V^2(K)\mathrm{Tr}\left(V^2(K)\right)E_d`$, where $`E_d`$ is the $`d`$-dimensional second-rank unit tensor) are of major interest for the physicist: the mass center of a rigid body moves in space as if the whole force acting on the body was concentrated onto a point particle at its center of mass. The inertia tensor determines how a body is rotating, given a fixed angular momentum. – Higher-rank moments such as symmetric $`r`$-rank tensors of the form $`V^r\left(K\right)dV𝐱^r`$ with $`r>2`$, may be considered, too.
Now it is again possible to write down a Steiner-type formula for the transition to the parallel body of a convex set $`K`$ beisbart\_schneider:schwerpunkte ; beisbart\_schneider:tensor : the scalar expansion coefficients have to be replaced by vector- or tensor-valued coefficients. The following simple argument shows that not only the character, but also the number of coefficients has to change: for a unit ball centered at the origin, $`B_R(0)\left\{𝐱\right||𝐱|R\}`$, the construction of the parallel body results in $`B_{(ϵ+R)}(0)`$. On the other hand, as a simple calculation shows, $`V^r\left(B_{(ϵ+R)}(𝐱)\right)\left(ϵ+R\right)^{d+r}`$. Thus, powers up to $`(d+r)`$ arise in the generalization of Steiner’s formula, which we write formally as

$$V^r\left(K_ϵ\right)=\underset{\nu =0}{\overset{d+r}{}}ϵ^\nu M_{r,\nu }.$$
(4)
As for the scalar Minkowski functionals, one can obtain an integral representation of the Steiner coefficients by considering a smooth and convex body $`K`$.

From the integral representation of the scalar MFs and from the form of the moments $`V^r`$ one might guess that the coefficients $`M_{r,\nu }`$ are built by functions of the form:

$$_K\mathrm{d}S^{d1}𝐱^rs_\nu (\kappa _1,..,\kappa _{d1})$$
(5)
with the elementary symmetric functions $`s_\nu `$ for $`\nu =0,..,(d1)`$ of the local curvatures $`\kappa _1,..,\kappa _{d1}`$ (for a definition, see beisbart\_schneider:tensor ). But this will not suffice, since we need $`(d+r)`$ coefficients. A closer look shows (see Figure 4), that indeed the local normal vector comes into play; for, if $`𝐱K`$, then $`𝐱+ϵ𝐧`$ lies on the surface of $`K_ϵ`$; here $`𝐧`$ denotes the local normal vector to $`K`$ at the point $`𝐱`$. The weighting with the local position in the definition of $`V^r`$ therefore leads to contributions from local normal vectors in the Steiner formula. In beisbart\_schneider:tensor  it is shown that in general the expansion coefficients $`M_{r,\nu }`$ can be written as linear combinations of the quantities:

$$W_\nu ^{p,q}\frac{1}{\nu \left(\genfrac{}{}{0pt}{}{d}{\nu }\right)}_K\mathrm{d}S^{d1}s_{\nu 1}(\kappa _1,..,\kappa _{d1})𝐱^p𝐧^q$$
(6)
for $`\nu =1,..,d`$, and $`p,q_0`$. Now the $`W_\nu ^{p,q}`$ look more basic than the coefficients arising in the Steiner formula for the volume tensor. It seems therefore, as if, in contrast to the case of the scalars, the Steiner coefficients are no longer the basic quantities. In order to generalize the Minkowski functionals we therefore have to move to an alternative route.

#### Spatial moments of the Minkowski functionals.

Let us take the moments $`V^r`$ together with the $`W_\nu ^{r,s}`$ ($`\nu =0,..,(d1)`$, $`r,s=0,1,..`$) as basic quantities. In order to arrive at a coherent notation, we set $`W_0^{r,0}V^r`$ and $`W_0^{r,j}=0`$ for $`j>0`$. Higher-rank Minkowski functionals therefore are higher moments of the scalar Minkowski functionals where the local surface is weighted with both position and surface normal. If the body $`K`$ is shifted with a translation vector $`𝐭`$, then the $`W_\nu ^{r,s}`$ transform like

$$W_\nu ^{r,s}\left(K+𝐭\right)=\underset{p=0}{\overset{r}{}}W_\nu ^{p,j}𝐭^{ip}\left(\genfrac{}{}{0pt}{}{i}{p}\right),$$
(7)
as follows immediately from their definition. For rotations $`R`$, we have

$$W_\nu ^{r,s}\left(RK\right)=RW_\nu ^{r,s}\left(K\right),$$
(8)
where we use the fact that the $`d`$-dimensional orthogonal group $`O(d)`$ has a natural representation in each tensor space.
Recalling that the Minkowski functionals can be described axiomatically, one would be inclined at this point to look for a characterization theorem for vectors and tensors. It turns out that one has to include even more valuations in order to achieve that purpose.

#### Axiomatic characterization.

A rigorous way to extend the Minkowski functionals is therefore to relax one of their defining axioms. We keep the additivity

$$\mathrm{\Phi }(AB)=\mathrm{\Phi }(A)+\mathrm{\Phi }(B)\mathrm{\Phi }(AB)$$
(9)
for a tensor valuation $`\mathrm{\Phi }`$; we also maintain the conditional continuity, i.e. continuity for the subclass of convex bodies222On the space of convex bodies one uses the Hausdorff metric, see beisbart\_schneider:brunn , p. 48 for a definition.; but, in the spirit of our previous findings, we replace the motion invariance by motion covariance defined as follows: for translations we require the existence of a number $`n_0`$ and of functions $`\mathrm{\Phi }_i`$ ($`i=0,..,n`$) such that for all $`K`$

$$\mathrm{\Phi }(K+𝐭)=\underset{i=0}{\overset{n}{}}\mathrm{\Phi }_i(K)𝐭^i.$$
(10)
Note, that products such as $`\mathrm{\Phi }_i(K)𝐭^i`$ always denote symmetric tensor products. For rotations we use the natural representation of $`O(d)`$ in the space of tensors as before, and require for $`RO(d)`$:

$$\mathrm{\Phi }(RK)=R\mathrm{\Phi }(K).$$
(11)
In a seminal work, Alesker beisbart\_alesker:tensor  proved the following characterization theorem: every motion-covariant, convex-continuous, and additive mapping $`\mathrm{\Phi }`$ can be written as a linear combination of the functionals

$$E_d^kW_\nu ^{r,s}$$
(12)
with $`i,j,k_0`$ and $`\nu =0,..,d`$. $`E_d`$ denotes the $`d`$-dimensional unit matrix333For higher than two-dimensional spaces, it does not make a difference whether we require $`O(d)`$\- or $`SO(d)`$-covariance. In the two-dimensional case, however, it does beisbart\_alesker:tensor . For the one-dimensional case see also beisbart\_alesker:rotation ..
Equation (12) gives all of the Minkowski valuations and finally answers the question how to generalize the MFs in a suitable manner. It contains more valuations than those given in Equation (6). However, some caution has to be exercised: Alesker’s theorem beisbart\_alesker:tensor  only bounds the dimension of the tensor functional spaces from above. It is not obvious that all of the above functionals are linear independent. Indeed, beisbart\_mcmullen:tensor  proved some linear relationships holding among the tensors. In order to give one example, Gauß’ theorem yields

$$W_1^{1,1}=_KdS^{d1}\mathrm{𝐱𝐧}=W_0^{0,0}(K)E_d.$$
(13)
One therefore has to explore the different tensor ranks further.

## 0.4 Exploring higher-rank Minkowski valuations

### 0.4.1 Some mathematical results

#### Minkowski vectors.

Minkowski vectors, i.e. Minkowski valuations of rank one, have been known for some years beisbart\_hadwiger:vekt ; beisbart\_schneider:schwerpunkte ; beisbart\_schneider:schwerpunkte2 ; beisbart\_hadwiger:vect2 . There is a strict one-to-one correspondence between the vectors and the scalars: for each Minkowski functional $`W_\nu _Ks_{\nu 1}(\kappa _1,..,\kappa _{d1})`$ there is a related vector $`_Ks_{\nu 1}(\kappa _1,..,\kappa _{d1})𝐱`$, such that $`(d+1)`$ vectors exist altogether. The reason why no normals enter at this level is basically because for compact convex bodies $`_Ks_{\nu 1}(\kappa _1,..,\kappa _{d1})𝐧=0`$ for $`\nu =0,..,(d1)`$ (for an elementary proof see beisbart\_hadwiger:vect2 ). The integral geometry of vectors was extensively scrutinized, see, e.g., beisbart\_hadwiger:vekt . For physical applications of vectors (see beisbart\_beisbart:querletter ; beisbart\_beisbart:fp ) the curvature centroids

$$𝐩_iW_\nu ^{1,0}/W_\nu ^{0,0}(\nu =0,..,d)$$
(14)
are a useful concept. Their meaning is straight forward: they locate morphological information. We shall see below that, if the curvature centroids do not coincide within one point, strong evidence for asymmetry is given.

#### Second-rank tensors.

The situation becomes more complicated for higher-rank tensors. It can be shown, that at the second rank, the following $`(3d+1)`$ tensors provide a base beisbart\_beisbart:tensor :

$`{\displaystyle _K}dV\mathrm{𝐱𝐱},`$ (15)
$`{\displaystyle _K}dS^{d1}s_\nu \mathrm{𝐱𝐱},{\displaystyle _K}dS^{d1}s_\nu \mathrm{𝐱𝐧},{\displaystyle _K}dS^{d1}s_\nu \mathrm{𝐧𝐧}`$ (16)
for $`\nu =0,..,(d1)`$. Not all of them are interesting from a physical point of view, since some of them are related to Minkowski functionals times the unit tensor and thus do not contain any new information. Principal kinematic formulae and Crofton formulae are known for second-rank tensors beisbart\_schneider:tensor ; beisbart\_beisbart:tensor ; beisbart\_schneider:tensor2 .
For a geometrical interpretation, we have to bear in mind, that the values of the tensors depend on the origin adopted. For $`W_\nu ^{2,0}`$ ($`\nu =0,..,d`$), however, the corresponding curvature centroid $`𝐩_\nu `$ is a natural reference point; the other tensors listed in Equation (16) are translation-invariant. Obviously, with such a choice the tensors are sensitive to morphological anisotropies, in the case of symmetry they align along the symmetry axes, otherwise they distinguish morphologically interesting directions. It is therefore useful to consider their eigenvalues and eigendirections.

#### Beyond second rank.

Although higher-rank tensors may be interesting as well, we shall not consider them further at this point.

### 0.4.2 Examples

In order to give a first example of how the Minkowski valuations discriminate spatial structures, we consider a few simple patterns as shown in Figure 5.

In the case of a circle (panel 1) the curvature centroids coincide within the circle’s center; the tensors – here represented as an ellipse featuring their eigenvalues and the corresponding directions – are isotropic. It is useful to divide the tensors $`W_\nu ^{r,s}`$ by the corresponding Minkowski functionals $`W_\nu ^{0,0}`$. The trace of the normalized tensor (i.e. the sum of the eigenvalues) quantifies to which extent the morphological information is concentrated around the center, how large the morphological fluctuations around the center are.
The transition to an ellipse of the same volume as the circle (panel 2) can be traced either with the tensors becoming anisotropic or with the ratio of the squared surface to the volume (isoperimetric ratio, see beisbart\_schmalzing:webI ). Furthermore, as the morphological concentration is decreasing, the traces of most of the normalized tensors increase, especially of $`W_2^{2,0}`$: the more elongated the ellipse becomes, the higher is its curvature at its very fringes along the horizontal direction. The normalized trace of $`W_1^{1,1}`$, on the contrary, becomes smaller, since the normal and the local position vector cease to be parallel to each other. The curvature centroids still coincide with each other because of the point symmetry. – If we reduce the symmetry further to a merely axially symmetric configuration (panel 3), the curvature centroids fan out along the symmetry axis. In the case of no symmetry at all (panel 4) they may constitute a triangle.
Finally, suppose one is given three convex bodies (or $`d+1`$ ones in $`d`$ dimensions) with known Minkowski functionals in order to combine them into a pattern such as indicated in panels 5 and 6. Given the scalar measures of the whole pattern, especially its Euler characteristic, one can constrain the number of intersections or subclumps in a pattern, which does not make any difference in our examples. For disentangling the positions of the grains, however, one has to know the curvature centroids. The centroids weight the subclusters according to their Minkowski functionals. Therefore, whereas from the volume’s point of view, only the larger components contribute significantly, all subclumps have the same weight (Euler characteristic $`=1`$) for the curvature centroid $`𝐩_2`$. If the subclumps do not overlap, one can indeed completely recover their positions.
In Figure 6, we show all tensors for the last example in more detail. Obviously, there are not only anisotropies; rather different morphological directions are relevant. The volume concentration around the center of mass is relatively high, since the volume is dominated by the two largest subclusters, being close to each other. It is quite the opposite with $`W_2^{2,0}`$: $`𝐩_2`$ is relatively far away from all the clumps, the curvature around $`𝐩_2`$ is dispersed leading to a high value of tensor’s trace. The main relevant direction connects the pair of large clumps with the third subcluster. The tensor $`W_1^{1,1}`$ is always isotropic, as shown in Equation (13), its trace is connected to the volume. Obviously, $`W_2^{1,1}`$ and $`W_1^{0,2}`$ complement each other; indeed they add up to a diagonal tensor (recall the linear relationships described in beisbart\_schneider:tensor ). Their eigenvectors’ alignment with the coordinate axes is due to the fact that all the subclumps are aligned parallel to the axes as well. Otherwise, $`W_2^{1,1}`$ and $`W_1^{0,2}`$ would average over the orientations of the subclumps. Finally, $`W_2^{0,2}`$ is isotropic, its trace equals the Euler characteristic.
The detailed description of patterns given by the Minkowski valuations can be used in practical applications, as we shall see in the next section.

## 0.5 Physical applications of higher-rank Minkowski valuations

### 0.5.1 The inner structure of spiral galaxies

The classification of spiral or disc galaxies into specific types (the Hubble types, e.g., beisbart\_hubble:classification ) is an important step towards their physical understanding; there are competing theoretical pictures trying to explain the formation of spiral arms beisbart\_lin:spiral ; beisbart\_gerola:stochastic . In Figure 7 we display a real galaxy observed face-on together with a simulation according to a percolation model beisbart\_gerola:stochastic . So far one has mainly relied on an assignment by eye in order to classify observed disc galaxies and to compare observational data with theoretical models beisbart\_naim:comparative . One of the reasons for this is that spiral arms are easily detected by the human eye, but are rather hard to quantify. Here we use Minkowski valuations to reveal aspects of the inner structure of spiral galaxies (see beisbart\_russell:spiral ; beisbart\_block:spiral ; beisbart\_kornreich:asymmetry  for some discussion of how to quantify morphometric features of spiral galaxies, and especially beisbart\_abraham:automated ; beisbart\_abraham:asymmetry ).
In order to render an optical galaxy image accessible to a Minkowski analysis, it is helpful to understand it as a two-dimensional surface brightness field $`u(𝐱)`$. We may then estimate the Minkowski valuations of the excursion sets $`M_v=\{𝐱|u(𝐱)>v\}`$. Contrary to large scale cosmological data (e.g. the distribution of galaxies or clusters of galaxies), such a density field is neither homogeneous nor isotropic; rather it reveals a point of considerable symmetry: the center of mass. We expect the curvature centroids to coincide within this point and therefore concentrate on second-rank Minkowski tensors.

In order to detect a possible overall anisotropy, we divide the smaller eigenvalue of $`W_0^{2,0}`$ by the larger one; this yields the parameter $`X`$. Since $`X1`$ by definition, $`ϵ1X`$ can serve as a measure of the anisotropy. Such an anisotropy may reflect a physical anisotropy, induced, for instance, by the presence of two massive spiral arms, but may also indicate that the galaxy is not observed face on, but rather under an inclination angle.
As the tensors $`W_\nu ^{2,0}`$ ($`\nu =0,..,2`$) weight the volume or the surface, respectively, with $`\mathrm{𝐱𝐱}`$, their eigenvalues are increased, if the configuration is more widespread (see Fig. 8). This fact can be used to quantify the concentration of a density field around its center of mass, by calculating the trace of $`W_\nu ^{2,0}`$.

Therefore another characteristic feature of disc galaxies, namely the way in which the radial surface brightness distribution declines (with increasing radius), can be analyzed by scanning through various threshold values $`v`$. We applied (suitably normalized) concentration parameters $`P_\nu `$ (built from $`\mathrm{Tr}(W_\nu ^{2,0})`$) to the images of both real galaxies (C. Moellenhoff, private communication) as well as to galaxies produced by a percolation model (beisbart\_seiden:percolation ; see beisbart\_dahlke:diploma  for details). These morphological descriptors reveal strong luminosity concentration for real galaxies, which could not be reproduced by any of the model galaxies (see Fig.9) with reasonable choices for the model parameters. Using the morphology of observed disc galaxies, one can therefore rule out the percolation model (at least in its simple form) as a viable attempt to explain the formation of spiral arms.

### 0.5.2 The morphology of galaxy clusters

As mentioned in Section 0.2, the morphology of clusters may constrain the values of the cosmological parameters. Intuitively, one may expect clusters in cosmological low-density models to have less substructure than in high-density models beisbart\_richstone:cluster . Using Minkowski valuations, we check whether this effect can be observed for the cosmological GIF simulations beisbart\_bartelmann:arcIV . We analyze the clusters’ projected mass distributions such as shown in Figure 2. Observations from the gravitational lensing through clusters beisbart\_schneider:book ; beisbart\_fischer:lens  trace approximately the same cluster information.
In order to make the simulation data, which consist of “particles tracing matter”, accessible to the Minkowski valuations, we construct a pixelized density field $`u(𝐱)`$ and smooth it with a Gaussian kernel of width444Length scales are given in units of $`h^1\mathrm{Mpc}`$. One Megaparsec (Mpc) equals about 3.26 million light years, the constant $`h0.65`$ accounts of the uncertainty in the measurements of the Hubble constant. $`\lambda =0.3h^1\mathrm{Mpc}`$. The resulting field $`u_\lambda (𝐱)`$ contains the overall shape of the cluster. Again we estimate the Minkowski valuations of the excursion sets, which contain all pixels above a given density threshold $`v`$. Our method to estimate the Minkowski functionals from the grid is based upon Crofton’s formulae beisbart\_serra:morphology ; beisbart\_schmalzing:beyond ; beisbart\_beisbart:diss . For each cluster we average a function of the Minkowski valuations $`f`$ over different density thresholds:

$$\frac{1}{\left(u_\lambda \right)_2\left(u_\lambda \right)_1}f_{\left(u_\lambda \right)_1}^{\left(u_\lambda \right)_2}du_\lambda f$$
(17)
in order to obtain a compact set of parameters. Here $`\left(u_\lambda \right)_1`$ and $`\left(u_\lambda \right)_2`$ denote appropriately chosen integration limits. We consider five types of substructure.

1. The clumpiness C $`\sqrt{\left(\frac{1}{\pi }W_21\right)^2}`$ counts the number of subclumps in quantifying how much the Euler characteristic $`W_2/\pi `$ deviates from one.
2. The construction of the asymmetry parameter A is based on the fact that curvature centroids which do not coincide with each other, indicate asymmetry. Therefore, we calculate the perimeter of the triangle constituted by the curvature centroids, $`2W_1\left(\mathrm{\Delta }\left(𝐩_i\right)\right)`$, A$`W_1\left(\mathrm{\Delta }\left(𝐩_i\right)\right)`$.
3. The ratio of the tensors’ eigenvalues, $`\gamma _1\left(W_\nu ^{r,s}\right)/\gamma _2\left(W_\nu ^{r,s}\right)`$ tells us to which extent anisotropies are present. Here we consider the tensor $`W_1^{2,0}`$ and define X$`\gamma _1\left(W_1^{2,0}\right)/\gamma _2\left(W_1^{2,0}\right)`$.
4. In general the curvature centroids $`𝐩_\nu `$ wander in space if the density threshold is varied; their fluctuations therefore measure the shift of morphological properties. Focusing on $`𝐩_1`$ we define S$`\left(𝐩_1𝐩_1\right)^2`$. This parameter is related to the centroid variation used in cluster investigations beforebeisbart\_mohr:x-ray-cluster .
5. Likewise, the eigendirections (parametrized by an angle $`\alpha `$) of the tensors are changing; this generates the twist of morphological properties. Restricting ourselves to $`W_1^{2,0}`$ we employ the parameter T$`\left(\alpha \alpha \right)^2`$.

For more technical details, see beisbart\_beisbart:diss .

For each time period (or cosmological redshift) we estimate the order parameters for three projections per cluster and average over all clusters in one cosmological model in order to extract a typical morphology. In Figure 10, we display the evolution of the averaged morphological parameters vs. redshift. One recognizes immediately a relaxation process, leading to less substructure (and lower values of the morphological order parameters) as time goes by (from right to left). Furthermore, clear differences between the cosmological background models arise. Indeed the theoretical predictions according to which the low-density clusters (OCDM) exhibit less substructure, are confirmed for the past; for redshift $`z=0`$, the differences disappear.
All of our order parameters show a similar behavior. The reason probably is that we average over different types of relaxation processes. An investigation of single clusters should reveal differences between, e.g., a bimodal merger dominated by two colliding subclumps and a more isotropic mass infall.

### 0.5.3 The geometry of the electric charge distribution in molecules

An important implication of the atom hypothesis is that the structures of molecules are to be explained in terms of their constituent parts, the atoms. Heuristic principles based upon this assumption work quite well in chemistry. However, if one tries to justify this hypothesis in the framework of quantum mechanics, the notion of an individual atom seems to disappear, since the negative charge distributions of the atoms merge, blurring the electron clouds which define the single atoms.
Bader beisbart\_bader:quantum  proposed an elegant method to recover the notion of atoms in molecules. Since, from the quantum mechanical point of view, the charge density profile $`\varrho (𝐫)`$ is the only relevant quantity, its geometry and topology must contain the notion of atoms. Using the critical points of the density distribution $`\varrho `$, where $`\varrho =0`$, one can define atoms as local density maxima together with a surrounding domain enclosed by zero-flux surfaces. This geometrical approach can be extended to chemical bonds and has been justified on physical grounds. Moreover, it provides us with means to visualize molecules in a proper way using the charge density profile, see Figure 11 for an illustration. A distinguished role is played by the surface $`^2\varrho =0`$ which separates the zone where chemical reactions are likely to happen from the rest of space.

Obviously, this is an area of research where physics and geometry strongly interact, and the applications of higher-rank Minkowski valuations may yield useful information on the structure of the chemical bonding. First, the morphology of the whole molecules can be simply described in terms of the Minkowski valuations of the density field. The atoms, which Bader defines beisbart\_bader:quantum , may not be compact in every case, so the single atoms can not be understood in terms of the Minkowski valuations. However, the fact that the atoms’ positions are specified in Bader’s theory gives us the chance of probing the molecular structure with the Boolean grain method where the atoms’ positions are decorated using spheres (either monodisperse or with radii mirroring physical quantities). Note, that at this point one really needs tensors, since in this case the locations and the orientations of subsystems are of crucial importance; and this local information lies beyond the scope of the scalar measures. Bader’s theory constructs also a graph of bonds, that may be analyzed with Minkowski valuations straightfowardly.
But given the fact that both the molecules and the constituting atoms are defined in terms of geometry, one may go one step further and try to discover phenomenological principles allowing one to predict the shape of molecules given their constituents.
These endeavors clearly are beyond the scope of this article. However, in order to illustrate to some extent how useful the Minkowski valuations may prove in this area, we consider briefly some of the Minkowski tensors for the H$`{}_{}{}^{+}{}_{2}{}^{}`$ ion.

In order to treat the H$`{}_{}{}^{+}{}_{2}{}^{}`$ ion consisting of two protons and an electron in the framework of quantum mechanics, one often considers the protons as fixed in space (Born-Oppenheimer approximation) thereby reducing the problem to an effective one-particle problem of an electron moving in the potential of two positively charged nuclei. A reasonable ansatz for the quantal wave function of the electron – at least suitable to illustrate what is basically happening – is a superposition of quantum states where the electron is bound to either of the protons. To be more precise, let us consider the protons as being separated by their equilibrium distance of approximately two Bohr radii555In the following we consider spatial distances in units of Bohr radii. The Bohr radius $`a_0`$ marks the typical scale of the hydrogen ground state wave function. For an introduction into the physical treatment of the H$`{}_{}{}^{+}{}_{2}{}^{}`$ ion one may consult, e.g. Complement GXI in beisbart\_cohen:qmii .. Let us label the protons with $`1`$ and $`2`$ and let $`\psi _i`$ denote the unperturbed (well-known) hydrogen ground state wave function of an electron bound to proton $`i`$ for $`i=1,2`$. For symmetry reasons, only the linear combinations

$$\psi _a=\frac{1}{\sqrt{2}}\left(\psi _1\psi _2\right)\mathrm{and}\psi _b=\frac{1}{\sqrt{2}}\left(\psi _1+\psi _2\right)$$
(18)
are of interest. Since $`\psi _b`$ is favored on energetic grounds, it is a bonding state; consequently $`\psi _a`$ is called antibonding state.
The morphology of the electron charge distribution is expected to distinguish between both states. Let us consider excursion sets of the three-dimensional electron charge density profile $`\varrho =|\psi |^2`$ for both states. In Figure 12 we display the surface tensor $`W_1^{2,0}`$ as a function of the charge fraction enclosed within the excursion set. We start from the very peaks of the charge distribution (high charge density thresholds) and move on to lower charge density thresholds enclosing a significant part of the whole electron charge. The trace of the tensor $`W_1^{2,0}`$ shows that the iso-surfaces of the bonding state are slightly more extended than those of the anti-bonding state. But on the whole, the differences are small: the concentration of the charge distribution is basically determined by the distance of the protons. If we ask for morphologically relevant directions, the differences between both states are much more perspicuous. The anisotropy of the azimuthal symmetric quantum states can be quantified via

$$An\left(\left(W_1^{2,0}\right)_{xx}\left(W_1^{2,0}\right)_{yy}\right)/\mathrm{Tr}\left(W_1^{2,0}\right),$$
(19)
the protons lying at $`(x,y,z)=(1,0,0)`$ and $`(1,0,0)`$, respectively. For small charge fractions, the tensor $`W_1^{2,0}`$ is dominated by the peaks around the positively charged protons yielding a maximum anisotropy, but no essential difference between both states. However, as we move to lower density thresholds and higher charge fractions, the bonding state isotropises much more quickly than the anti-bonding one. In the bonding state, therefore, the electron forms a cloud spreading out over both atoms in a quite isotropic way, whereas in the anti-bonding state, the signature of the two distinct protons remains still visible. The chemical bond therefore is mirrored by the geometry of the electron density distribution, as quantified via Minkowski tensors.
This example is but an elementary first illustration of the nature of the chemical bond. In the spirit of Bader’s theory, one can possibly go much further in order to consider single atoms in molecules. Additivity can serve as a guiding principle, since Bader found that there are stable atom configurations which persist through different types of molecules, presenting a sort of rigid module movable as a whole; this fact is well-known from chemistry, of course.

## 0.6 Density functional theory based on Minkowski valuations

Having shown the morphometric versatility of MVs, we conclude our discussion of applications by considering systems where integral geometric quantities enter a physical approximation. The quantities that come into play at this point generalize the Minkowski valuations further and may indicate a new direction of interest.
Density functional theory became the standard method to study inhomogeneous classical fluids in the last 20 years. The main feature of such physical systems is the spatial variation of the equilibrium density $`\rho _{eq}(𝐱)`$ which minimizes the free energy $`\mathrm{\Omega }`$. Density functional theory is based on the exact result that the free energy of the inhomogeneous fluid can be expressed as a functional $`\mathrm{\Omega }[\rho (𝐱)]`$ of the density $`\rho (𝐱)`$. Thus, all relevant thermodynamic properties such as surface tensions or phase transitions of confined fluids can be calculated from $`\mathrm{\Omega }[\rho _{eq}(𝐱)]`$; even equilibrium correlation functions that describe the microscopic structure of inhomogeneous fluids can be determined by derivatives of the functional.
In general, the form of the functional $`\mathrm{\Omega }[\rho (𝐱)]`$ is not known; deriving the exact functional would be equivalent to solving the statistical mechanics for fluid systems exactly which seems to be an impossible task. Therefore, one has to restrict oneself to reasonable approximations for $`\mathrm{\Omega }[\rho (𝐱)]`$. For convenience, one separates an excess intrinsic free energy functional $`\mathrm{f}_{ex}[\{\rho _i\}]`$ from exactly known external contributions

$$\mathrm{\Omega }[\rho (𝐱)]=d𝐱\left(\mathrm{f}_{ex}[\rho (𝐱)]+k_BT\rho (𝐱)\left(\mathrm{ln}\mathrm{\Lambda }^3\rho (𝐱)1\right)+\rho (𝐱)\left(V(𝐱)\mu \right)\right)$$
(20)
where $`V(𝐱)`$ is an external potential, $`\mu `$ the chemical potential, $`\beta =1/(k_BT)`$ the inverse temperature $`T`$, and $`\mathrm{\Lambda }`$ the thermal wavelength of the fluid. $`k_B`$ denotes the Boltzmann constant.
The most elaborated and reliable functional for hard-sphere fluids is that developed by Rosenfeld beisbart\_rosenfeld:free  based on the local scalar and vector-valued Minkowski functionals

$$W_\nu ^{p,q}[\rho (𝐱)]\frac{1}{\nu \left(\genfrac{}{}{0pt}{}{d}{\nu }\right)}_K\mathrm{d}S^{d1}s_{\nu 1}(\kappa _1,..,\kappa _{d1})𝐱^p𝐧^q\rho (𝐱).$$
(21)
In contrast to the definition (6) and the applications considered in the previous sections, the functionals here are locally weighted with the density $`\rho (𝐱)`$. Notice, that for inhomogeneous densities the vector valuations $`W_\nu ^{0,1}[\rho (𝐱)]`$ in general do not vanish for a convex body $`K`$ and that the functionals do inherently depend on the location of $`K`$. The following ansatz for the excess free energy functional proves to be useful:

$$\begin{array}{c}\beta \mathrm{f}_{ex}[\rho (𝐱)]=f_1(W_0^{0,0})W_3^{0,0}+f_2(W_0^{0,0})W_1^{0,0}W_2^{0,0}+f_3(W_0^{0,0})(W_2^{0,0})^3\\ \\ +f_4(W_0^{0,0})W_1^{0,1}W_2^{0,1}+f_5(W_0^{0,0})W_2^{0,0}(W_2^{0,1}W_2^{0,1}),\end{array}$$
(22)
where we drop the argument $`[\rho ]`$ for convenience. This ansatz is inspired by the observation that the physical dimensions are $`[W_3^{0,0}]=[W_1^{0,q}W_2^{0,q}]=[(W_2^{0,q})^3]=[\beta f_{ex}]=[\mathrm{length}]^3`$, that $`[W_0^{0,0}]=1`$ is a dimensionless weighted density.
To obtain the functions $`f_i(W_0^{0,0})`$ one may apply an argument from scaled particle theory for a homogeneous fluid beisbart\_reiss:rigid , namely, that the excess chemical potential for big spheres of radius $`\sigma /2\mathrm{}`$ is $`lim\mu /Vp`$ with the volume $`V=\frac{\pi }{6}\sigma ^3`$, the chemical potential $`\mu =\frac{\mathrm{f}_{ex}}{\rho }`$ and the pressure $`p`$. Within the present ansatz (22) one finds in this limit the reversible work to create a cavity of volume $`V`$, $`\mu =_\nu \frac{\mathrm{f}_{ex}}{W_\nu ^{0,0}}\frac{W_\nu ^{0,0}}{\rho }`$, and therefore in the limit $`\sigma /2\mathrm{}`$ the chemical potential $`\mu /V\frac{\mathrm{f}_{ex}}{W_0^{0,0}}`$. Using the equation of state $`\mathrm{\Omega }=pV`$, i.e., $`p=\mathrm{f}_{ex}+_\nu \frac{\mathrm{f}_{ex}}{W_\nu }W_\nu +W_3^{0,0}`$ one gets the well known ’scaled particle’ differential equation

$$\frac{\mathrm{f}_{ex}}{W_0^{0,0}}=\mathrm{f}_{ex}+\underset{\nu =0}{\overset{3}{}}\frac{\mathrm{f}_{ex}}{W_\nu ^{0,0}}W_\nu ^{0,0}+\frac{\mathrm{f}_{ex}}{W_1^{0,1}}W_1^{0,1}+\frac{\mathrm{f}_{ex}}{W_2^{0,1}}W_2^{0,1}+W_3^{0,0}.$$
(23)
Within the ansatz (22) the solution of the differential equation reads

$$\begin{array}{cc}\beta f_{ex}[\rho (𝐱)]=\hfill & \frac{3W_3^{0,0}}{4\pi }\mathrm{ln}(1W_0^{0,0})+\frac{9}{4\pi }\frac{W_1^{0,0}W_2^{0,0}W_1^{0,1}W_2^{0,1}}{1W_0^{0,0}}\hfill \\ & \\ & +\frac{9W_1^{0,0}}{8\pi }\frac{(W_1^{0,0})^23W_1^{0,1}W_1^{0,1}}{(1W_0^{0,0})^2}.\hfill \end{array}$$
(24)
This functional yields excellent results for thermodynamic quantities and for structures in the fluid phase of hard spherical molecules beisbart\_roth:depletion . For instance, in Figure 13 the inhomogeneous density profile $`\rho (z)`$ close to a hard wall is shown as compared to simulation results.
Unfortunately, the functional in Eq. (24) was derived only for hard spheres in three

dimensions; however, a proposal has been made for general convex bodies beisbart\_rosenfeld:molecular  based on an approximative decomposition of $`f_{ex}`$ in the low density limit. Here, we show how Rosenfeld’s idea of decomposing $`f_{ex}`$ for spheres beisbart\_rosenfeld:free  may be extended exactly to non-spherical convex bodies by applying integral geometric techniques such as kinematic formulae and an explicit local expression for the Euler characteristic $`\chi (KK^{})`$ of two overlapping bodies beisbart\_mecke:fundamental . The main idea is to decompose the exact expression for the local excess free energy $`f_{ex}[\rho (𝐱)]`$ at low densities in terms of the vector-valued Minkowski functionals. To this end one has to evaluate the configurational integral $`dK=dRd𝐱`$ over all orientations $`R`$ and positions $`𝐱`$ of a body $`K=RK_𝐱`$, namely,

$$\begin{array}{cc}& \frac{4\pi }{3}d𝐱\beta f_{ex}[\rho (𝐱)]\frac{2\pi }{3}dKdK^{}\chi (K_𝐱K_𝐱^{})\rho (𝐱)\rho (𝐱^{})\hfill \\ & \\ =& _^dd𝐱W_3^{0,0}[\rho (𝐱)]W_0^{0,0}[\rho (𝐱)]\hfill \\ +& _^dd𝐱\mathrm{\hspace{0.33em}3}W_1^{0,0}[\rho (𝐱)]W_2^{0,0}[\rho (𝐱)]3W_1^{0,1}[\rho (𝐱)]W_2^{0,1}[\rho (𝐱)]\hfill \\ +& _^dd𝐱_𝐫dK_𝐱\rho (𝐫)_𝐫^{}dK_𝐱^{}\rho (𝐫^{})\frac{\kappa _1\kappa _2}{2}\frac{\left(𝐯_1𝐯_1𝐯_2𝐯_2\right)_{ij}(𝐧^{}𝐧^{})_{ij}}{1+\mathrm{𝐧𝐧}^{}}\hfill \end{array}$$
(25)
where $`𝐯_i`$ denotes the principal direction of the curvature $`\kappa _i`$ at $`𝐫`$ and $`𝐧`$ ($`𝐧^{}`$) are the normal vectors of $`RK_𝐱`$ ($`R^{}K_𝐱`$) at $`𝐫`$ ($`𝐫^{}`$), respectively. The integral over rotations $`R`$ and over the surface of the body $`K_𝐱`$ centered at $`𝐱`$ is written as $`_𝐫dK_𝐱=dR_{𝐫RK_𝐱}dS^{d1}`$. For spheres with $`\kappa _1=\kappa _2`$ one recovers immediately the low-density limit of the intrinsic free energy given by Eq. (24). But unfortunately the last term on the right hand side of Eq. (25) vanishes only for spheres, so that only for hard spheres a finite number of Minkowski functionals is, by accident, sufficient to decompose the excess free energy into contributions stemming alone from the individual bodies, whereas for non-spherical shapes in arbitrary dimensions an infinite number of fundamental curvature measures is required. Nevertheless, based on the most general analytical expression (25), several approximations by a finite number of tensor-valued Minkowski functionals (including the proposal in beisbart\_rosenfeld:molecular ) may be proposed and numerically tested which make Rosenfeld’s density functional applicable to general convex bodies, in particular, to spherocylinders and ellipsoids beisbart\_mecke:fundamental .

## 0.7 Conclusions

Concluding, we argue that a new perspective emerges as regards the Minkowski functionals. Instead of marking a small number of distinguished measures (which, of course, they do in a certain sense), the Minkowski functionals are embedded into a more extensive class of Minkowski valuations. These measures generalize the Minkowski functionals without leaving the framework of integral geometry. The essential step to extend the scalar functionals is to replace the motion invariance by motion covariance. This step yields a hierarchy of tensor-valued measures with the Minkowski functionals at the bottom. The higher-rank valuations are moments of the Minkowski functionals where local geometry is weighted with either the position or the normal vector. The full, comprehensive description of the tensor-valued measures is still a mathematical challenge, although some important results concerning the integral geometry of tensor measures have already been achieved, especially in low dimensions (see beisbart\_schneider:tensor2 ).
The extensions of the Minkowski functionals are not only mathematically interesting; rather they can be motivated on physical grounds. Since the Minkowski functionals proved to be useful for a number of applications, we are justified to expect that also their generalizations may be suitable tools. Using physical applications at various length scales we exemplified how higher-rank Minkowski valuations complement the information featured by the scalar functionals. From this point of view we hope that the higher-rank Minkowski valuations will become a standard tool for morphometry and some other fields of research concerned with pattern formation.

## Acknowledgement

We thank C. Möllenhoff for providing us with recent observations of spiral galaxies. Furthermore, we thank C. Levit for the permission to show the images of the molecules and J. Colberg for the GIF cluster data. We thankfully acknowledge comments on the manuscript by M. Kerscher and R. Schneider. Last but not least we thank the participants of the second Wuppertal conference for vivid discussions.


format: text

text:
# Topological (Sliced) Doping of a 3D Peierls System: Predicted Structure of Doped BaBiO3

## I introduction

When a 1D crystal is driven incommensurate, for example, by doping to alter the Fermi wavevector, it is well-known that the modulation of the commensurate crystalline order tends not to be uniform but to accumulate at kinks in the crystalline order parameter . A prototype for this is the Su-Schrieffer-Heeger model  for polyacetylene, (CH)x. In a simple chemical view, the carbon chain has alternating single and double bonds. Crystallography confirms this lattice dimerization. The scalar order parameter (amplitude of the staggered charge disproportionation) has two possible values, $`\pm \rho `$. Zero-dimensional domain walls (topological defects) separate regions of the 1D chain with positive and negative order parameters. Su et al. showed that when carriers are doped into polyacetylene, they create new domain-wall defects and localize into mid-gap electronic soliton states on these defects. Analogous effects in two dimensions  occur in layered materials near a Mott insulator phase, such as cuprates, nickelates, and layered manganites. The one-dimensional domain wall in the 2D system is generally called a “stripe.” The 3D analog, which occurs in some cubic manganites , is a planar domain wall , sometimes called a “sheet” or “lamella.” We prefer the term “slice,” which, like the word “stripe” can be used as verb or noun. Here we discuss a 3D model system in which bonding and strain energies favor slices, but long-range Coulomb interactions prefer point defects. We analyze the system numerically and show that at doping levels near 30%, Ba1-xKxBiO3 is near the crossover where slices may become more favorable than point defects for lowering the energy of the doped-in holes.

Charge-density-wave (CDW) and Peierls systems have scalar order parameters $`\rho _Q=N^1_{\mathrm{}}\rho _{\mathrm{}}\mathrm{exp}(i\stackrel{}{Q}\stackrel{}{\mathrm{}})`$. BaBiO3 is a simple 3D example. With nominal valence Bi4+, the Bi 6$`s`$-band is half-filled. The material is non-metallic, with a 2eV optical gap . Crystallography  shows a doubled unit cell; $`\rho _{\mathrm{}}`$ alternates with wavevector $`Q=(\pi ,\pi ,\pi )`$. The simple cubic sublattice of nominal Bi4+ ions self-organizes into a bipartite (rocksalt-type) charge-ordered array of nominal Bi3+ ($`\rho =4\rho _Q`$, called “A” sites) and nominal Bi5+ ($`\rho =4+\rho _Q`$, called “B” sites) ions. The actual value of the order parameter $`\rho _Q`$ has magnitude $`1`$, and takes two degenerate values, $`\pm |\rho _Q|`$. A simple way to think of this is as a 3D version of a Peierls instability. In a bi-partite lattice with only nearest-neighbor hopping there is accidental Fermi surface nesting: $`ϵ(\stackrel{}{k})=ϵ(\stackrel{}{k}+(\pi ,\pi ,\pi ))`$, and both are zero at the Fermi energy. This guarantees that the crystal can reduce its electronic band energy via the electron-phonon interaction by dimerizing. A Hamiltonian which contains this effect was introduced by Rice and Sneddon .

Regions of charge-ordered BaBiO3 with $`\rho _Q>0`$ are separated from regions with $`\rho _Q<0`$ by stacking faults. The simplest stacking fault lies in a (111) plane. In perfectly ordered BaBiO3, (111) planes are alternating A and B type (Fig. 1a). A stacking fault with no nuclear disorder has either two adjacent A layers (local charge excess -1 per site on the plane, or electron-doped) or two adjacent B layers (local charge excess 1, or hole-doped) as in Fig. 1b.

Here we point out that the Rice-Sneddon Hamiltonian, given below, predicts that holes or electrons, when introduced, will self-organize into slices. However, the long-range Coulomb interaction, neglected in the Rice-Sneddon model, will destabilize slices except possibly at the most favorable doping level. For light doping, the preferred structure for holes is to self-organize into point defects, small bipolarons, as previously discussed .

## II Rice-Sneddon Model

The Rice-Sneddon model  is simple, well-studied , and quite successful . In the perovskite crystal structure of BaBiO3, each Bi atom is surrounded by six oxygens, and each oxygen is shared by two Bi atoms. At high temperatures (of order 1000K) the crystal is nearly cubic perovskite, but at lower $`T`$ there are rotations and distortions of the BiO6 octahedra, which enlarge the unit cell. We believe that the most important effect is the “breathing” displacements $`u_0\pm 0.12\AA `$  of oxygens away from or toward alternate Bi atoms. These provide a natural electron-phonon mechanism to enhance the Bi 6$`s`$ charge density on A sites where oxygens breath outward, and reduce the Bi 6$`s`$ charge density on B sites where the oxygens breathe in. A microscopic Hamiltonian containing the minimal necessary electron-phonon interaction was given by Rice and Sneddon ,

$`=`$ $``$ $`t{\displaystyle \underset{<\mathrm{},\mathrm{}^{}>}{}}c_{\mathrm{}}^{}c_{\mathrm{}^{}}g{\displaystyle \underset{\mathrm{}}{}}e(\stackrel{}{\mathrm{}})c_{\mathrm{}}^{}c_{\mathrm{}}`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}K{\displaystyle \underset{\mathrm{},\alpha }{}}u(\stackrel{}{\mathrm{}}+{\displaystyle \frac{\widehat{\alpha }}{2}})^2.`$ (2)
The first term is nearest-neighbor hopping of Bi $`6s`$ electrons with hopping integral $`t`$ 0.35 eV. The index of summation $`\mathrm{}`$ implicitly includes a spin as well as site quantum number. The filling is $`1x`$ electrons per site. The variable $`u(\stackrel{}{\mathrm{}}+\widehat{\alpha }/2)`$ (with $`\alpha =x,y,z`$) is the displacement along a Bi-O-Bi bond in the $`\widehat{\alpha }`$ direction of the oxygen located at position $`(\stackrel{}{\mathrm{}}+\widehat{\alpha }/2)a`$. The variable $`e(\mathrm{})`$ is the local dilation or “breathing” amplitude of the 6 oxygens which surround the Bi ion at site $`\mathrm{}`$,

$`e(\mathrm{})`$ $`=`$ $`[u(\stackrel{}{\mathrm{}},x+\widehat{x}/2)u(\stackrel{}{\mathrm{}}\widehat{x}/2)]+[xy]`$ (3)
$`+`$ $`[xz].`$ (4)
The Einstein restoring force $`K`$ 19 eV/$`\AA ^2`$ is fitted to the measured 70 meV frequency of the Raman-active Peierls breathing mode .

At half-filling, this model opens a Peierls gap at the Fermi level; the electron-phonon interaction parameter $`g`$ 1.39 eV/$`\AA `$ is fitted to the measured  $``$2 eV gap. This is an ordinary size electron-phonon coupling. The change in Coulomb field of a charge -2$`e`$ oxygen ion gives $`g`$ 1.2 eV/$`\AA `$. The resulting dimensionless coupling constant $`\mathrm{\Gamma }g^2/Kt`$ is $``$0.30, intermediate between the weak ($`\mathrm{\Gamma }<0.2`$) and strong ($`\mathrm{\Gamma }>0.4`$) coupling regimes. In this middle regime, neither hopping nor electron-phonon energy dominates . The ground-state of undoped BaBiO3 is as close to a bipolaronic crystal (large $`\mathrm{\Gamma }`$, $`|\rho _Q|1`$) as to the conventional Peierls-CDW (small $`\mathrm{\Gamma }`$). We calculated  the order parameter $`\rho _Q`$ at $`\mathrm{\Gamma }=0.3`$ to be 0.82. The corresponding oxygen displacement $`u_0=2g\rho _Q/K=0.12\AA `$ agrees with the diffraction and EXAFS measurements , showing that the model is internally consistent.

We previously found that excitations across the Peierls gap form self-trapped excitons . We also reported that holes inserted into BaBiO3 self-trap and form bipolarons  since the coupling strength exceeds $`\mathrm{\Gamma }_c=0.17`$. These are doubly charged point defects, corresponding to local depressions of the order parameter where the oxygen distortion $`e(\stackrel{}{\mathrm{}})0`$ for $`t=0`$. For non-zero hopping $`t`$ the bipolaron spreads out and evolves continuously from a small bipolaron ($`\mathrm{\Gamma }\mathrm{\Gamma }_c`$) to large CDW-like bipolaron as $`\mathrm{\Gamma }\mathrm{\Gamma }_c`$. The stability of bipolaron defects provides a simple explanation why dilutely doped BaBiO3 remains insulating and diamagnetic.

The “disproportionation reaction” 2Bi$`{}_{}{}^{4+}`$Bi3++Bi5+ has been much discussed in the literature on BaBiO3. Using the definition $`2U=E(`$Bi$`{}_{}{}^{3+})+E(`$Bi$`{}_{}{}^{5+})2E(`$Bi$`{}_{}{}^{4+})`$, one can say that the effective Hubbard $`U`$ parameter is negative. Two factors may contribute to $`U`$, electron-phonon effects (expected to be attractive) and Coulomb repulsion ($`U=U_{\mathrm{ep}}+U_{\mathrm{el}}`$). If one wants to assign the mechanism for disproportionation completely to Coulomb interactions, then the Hubbard $`U`$ calculated with all atoms held stationary in cubic perovskite positions (defined as $`U_{\mathrm{el}}`$) should be negative. Vielsack and Weber  did careful calculations of $`U_{\mathrm{el}}`$, finding no evidence for negative values, but instead a small positive value $`U_{\mathrm{el}}0.6\pm 0.4`$eV. Therefore we shall temporarily ignore the on-site Coulomb repulsion $`U_{\mathrm{el}}`$. Apparently the Bi4+ ion must reorganize its environment in order to stabilize the charge disproportionation.

## III Hole doping

What happens at finite doping concentrations $`x`$? Assuming sufficient electron-phonon coupling to destroy the undistorted metal, there are two possibilities, (1) bipolarons, and (2) slices.

(1) Numerical studies by Yu et al. , and confirmed by us, show that randomly located point bipolarons are at least metastable. Depending on whether bipolarons attract or repel, the system could then either phase-separate into undoped and doped regions, or form spatially separated bipolarons. When bipolarons are small, the energy of an array of bipolarons is described by a pair-wise additive potential $`V(\mathrm{\Gamma },r)`$, containing the repulsive long-range Coulomb interaction $`V_{\mathrm{Coul}}`$ (neglected for the time being) and the interaction $`V_t`$, induced by hopping. $`V_t`$ decays exponentially with the distance $`r`$ between two bipolarons (like bipolaron wave functions.) Since bipolarons can only sit on former A sites, the nearest neighbor interaction $`V_0=V_t(\mathrm{\Gamma },\sqrt{2}a)`$, which is also the strongest interaction, is between bipolarons separated by $`a(110)`$.) We computed $`V_0`$ numerically by optimizing the oxygen-positions self-consistently for given bipolaron positions. Fig. 2 shows that large CDW-like bipolarons attract each other, whereas small bipolarons repel, with $`V_t0`$ in the atomic limit $`t0`$. Perturbation theory around $`t=0`$ gives a bipolaron repulsion. At the physically relevant value $`\mathrm{\Gamma }=0.3`$, bipolarons repel according to Fig. 2, but multi-bipolaron interactions become important with decreasing $`\mathrm{\Gamma }`$. Therefore, in the intermediate coupling regime we must rely on exact numerical diagonalization of Eq. (2) .

(2) In contrast to bipolaron defects where the order parameter never changes sign, holes could form topological defects. Phase-slips are planar defects with sign-changes of the charge order parameter $`\rho _Q`$ and of the breathing order parameter $`\widehat{e}=(1)^le_l`$. Consider a “BB” stacking fault in the 111-direction (Fig. 1b). In the atomic ($`t0`$) limit, each B-site on a phase-slip plane has three displaced and three undisplaced oxygen neighbors, i.e. $`e_l\pm 3u_0`$. Phase-slips accumulate one hole for every two atoms on a (111) BB bilayer. The average hole charge on phase-slip B-sites is thus $`\rho _{\mathrm{hole}}+\frac{1}{2}`$. The actual values of hole charge found for $`\mathrm{\Gamma }=0.3`$ and $`x=1/4`$ are shown in Fig. 3.

Our aim is to find the ground-state hole arrangement, testing the stability of several bipolaronic (1) versus the phase-slip (2) solutions numerically. We did a series of calculations on $`x=1/4`$-doped systems. The bipolaron sytems (1) were (a) maximal spacing between bipolarons, obtained when they occupy center and corner sites of a tetragonal 16 atom unit cell (bct structure); (b) a simple cubic (sc) arrangement of bipolarons sitting at the corner sites of a 8 atom cubic cell; (c) a disordered structure with random bipolaron positions; (d) phase-separated structures based on unit cells containing 8 or 16 (111) planes, where we replace one or two near neighbor A planes with bipolaronic B planes, \[(AB)3B2\]N or \[(AB)6B4\]N/2. We also looked at unit cells with 8 or 16 111 planes containing phase slips (2) \[AB(ABB)2\]N and \[(AB)2(ABB)4\]N/2. Finally, we looked at the undistorted metal \[C\]N where each atom C has a nominal Bi4.25+ valence.

Phase separation was strongly disfavored, while separated bipolarons and phase slips were all metastable. The phase slips weakly repelled, preferring the 8 plane solution to the 16 plane solution. The order parameter and hole charge density of this 8 plane solution are shown in Fig.(3). The stability is determined computing the total energy $`E_{\mathrm{tot}}(\{u_\stackrel{}{l}\})`$ given by Eq.(2), which is a function of the oxygen-positions $`u_\stackrel{}{l}`$. We start by guessing oxygen positions to get close to a local minimum in the energy landscape, then vary oxygen positions using a gradient minimization routine to find a self-consisistent minimum. For smaller periodic structure, we used $`k`$-space sampling in the corresponding Brillouin zones (8000 k-points). For each $`\stackrel{}{k}`$, $``$ is diagonalized exactly. States are filled with two electrons up to the desired doping. For the random bipolaron structure we used large asymmetric clusters ($``$400 atoms) with periodic boundary conditions and $`k=0`$ only. Initial oxygen-positions had Peierls order with small random deviations. Our calculations on random configurations generally reproduce the earlier calculations of Yu et al. .

Energies of various states at doping $`x`$=1/4 are shown in Fig. 4. In the atomic limit $`t=0`$ all bipolaron and phase-slip structures are degenerate. All that is required is that no A site should have an A first neighbor. Below a critical coupling strength $`\mathrm{\Gamma }_c(x)`$ ($`0.2`$ for doping $`x`$=1/4), distorted structures become unstable with respect to the undistorted metallic structure. Above $`\mathrm{\Gamma }_c(x)`$, we find numerically that at $`x=1/4`$ holes strongly prefer to order in phase-slips (2). Bipolaronic structures (1) are quite similar to each other in energy, and behave as expected from the the bipolaron-pair interaction. In the intermediate coupling range there occurs a cross-over from the tetragonal to the layered BBB-structure as stable bipolaron configuration, corresponding to a change in the overall bipolaron interaction from weakly repulsive to weakly attractive. For 1/4-doped BaBiO3, the energy gain for the phase-slip solution compared to the bipolarons is about 50-65 meV per hole.

Thus, contrary to previous studies of doping of this model , we find that the stable doping state is not bipolarons but phase slips. Does this model correspond sufficiently to reality for BaBiO3? Iwano and Nasu  use a more complicated hopping and electron-phonon interaction, but we believe that such corrections are not the relevant ones. There are also (i) small structural distortions (rotations of oxygen octahedra) beyond the breathing-mode distortions considered here, and (ii) non-adiabatic effects (such as zero point motion) associated with the fact that the oxygen mass is not infinitely large compared with the electron mass. We believe that both of these also have little relevance. It is harder to dismiss two other effects: (iii) the disorder caused by the dopant atoms, and (iv) the long-range Coulomb interaction, both omitted so far. Of these, the last is clearly important, as we now show, and tends to destabilize phase-slips relative to distributed bipolarons.

## IV Long-range Coulomb effect

At low doping, there is a large Coulomb cost in putting charges onto stacking faults instead of widely distributed point charges. We modeled this as follows. The Madelung energy was computed by the Ewald method for (Ba2+)2(Bi3+Bi5+)(O2-)6. The calculation was re-done for many large unit cells, with holes added on selected Bi ions (Bi$`{}_{}{}^{3+}`$Bi5+), and compensating negative charges distributed uniformly throughout space. This is one way to mimic potassium doping of BaBiO3. First consider the sliced state. By numerical calculation for uniformly distributed phase slips at many values of $`x`$ between 1/39 and 1/3, we found a good fit to the formula $`E_S=(25.1+\pi /9x)(e^2/2aϵ_{\mathrm{}})`$. The static electronic screening $`ϵ_{\mathrm{}}5`$ was measured by Tajima et al. . $`E_S`$ is the difference of energy per hole between the sliced solution and the undoped Peierls insulator. The term $`\pi /9x`$ is the analytic result for idealized uniform sheets of charge of vanishing thickness, arranged periodically in a compensating charge background. The term -25.1 corrects for the discreteness of the charges, the absence of the self-interaction, and includes the Coulomb energy of the holes with the background BaBiO3 lattice.

We also need the energy difference per hole $`E_B`$ of optimally spaced bipolarons relative to the undoped Peierls insulator. For doping $`x=1/n^3`$ the bipolarons can be placed on a sublattice of face-centered cubic (fcc) form, with maximum spacing. Numerical results for $`n`$=2, 3, and 4 fitted well to the formula $`E_B=(21.24.585x^{1/3})(e^2/2aϵ_{\mathrm{}})`$. The term $`4.585x^{1/3}`$ is the Madelung energy of an fcc lattice of charges $`2e`$ in a uniform compensating background, and the constant -21.2 accounts for the energy of the holes in the background BaBiO3 lattice. The difference $`E_SE_B`$ is the Coulomb penalty per hole for forming slices or phase slips. It is plotted in Fig 5, becoming large at low $`x`$ with a shallow minimum near $`x=0.3`$.

A more accurate estimate of the Coulomb penalty is shown as an x on Fig 5, and was obtained by doing Ewald sums using the actual Bi site charges ($`\rho _i=5n_i`$ where $`n_i=2<\psi _i^{}\psi _i>`$ is the local occupancy of the Bi $`s`$ states) and the actual (ABBABABB)N slice structure and bct bipolaron structure. It is seen on the figure that at 1/4 doping, the Coulomb penalty is twice as high as the Peierls benefit in forming slices. These calculations did not account for the actual random positions of the compensating negative charges where Ba2+ ions are replaced by K+ ions. Instead, the compensating charge was distributed uniformly. The calculations were repeated for a model where the compensating negative charges were equally shared by each Ba atom (Ba$`{}_{}{}^{2+}`$Ba(2-x)+). This only changed the constant terms in $`E_S`$ and $`E_B`$, but did not affect the difference, $`E_SE_B`$ plotted here.

The Madelung sums discussed above are not the complete Coulomb effect. The missing on-site repulsion has the opposite effect, preferring sliced solutions to distributed bipolarons. If we add back the local term

$$_U=U_{\mathrm{el}}\underset{i}{}(\rho _i4)^2$$
(5)
and treat it as a first-order perturbation, the sliced solution at 1/4-doping has lower on-site energy by 0.121$`U_{\mathrm{el}}`$ per hole than the bct bipolaron solution. Using the value $`U_{\mathrm{el}}0.6\pm 0.4`$eV from Vielsack and Weber , the onsite correction $`0.07\pm 0.05`$eV is potentially sufficient to re-stabilize the sliced solution. Thus we can say that the best place to look for sliced structures in BaBiO3 is in the range near 1/4 to 1/3 doping, but we cannot predict whether the sliced solution will be destroyed by disorder or Coulomb effects.

Doping 1/3 is the highest at which a simple sliced solution ($`(ABB)_N`$) is possible. For higher doping, either the A-planes acquire polaron defects or the sliced solution is destroyed and purely bipolaronic states win.

It is worth mentioning that there have been reports  (subsequently attributed to electron beam heating effects ) of superlattice diffraction spots in doped BaBiO3, not apparently identical to the superstructures predicted here. A further search in the doping region $`0.2<x<0.35`$ would be interesting.